Analysis of a Vector Field and Some 
Applications to Fluid Motion 


W. J. DUNCAN 


(University of Glasgow) 


SUMMaRY: The paper first discusses the analysis of a general vector field with 

the aid of the concept of the field line with its tangent, principal normal and 

binormal. Detailed applications are made to the motion of inviscid and of 
viscous fluids. 


1; Introduction 


The first part of this paper gives a general analysis of a vector field in which 
the basic concept is the field line; this covers the important special cases of line of 
force, streamline and vortex line. Simple relations are obtained for the rates of 
change of certain quantities along the tangent, principal normal and binormal of 
the field line. 


The applications to fluid motion are based on identifying the general vector of 
the theory with the instantaneous velocity vector or with the instantaneous vorticity 
vector. The applications are most easy and valuable when the motion is steady 
and cover an extension of Bernoulli’s theorem to viscous fluids. 


It is not claimed that the conclusions of the paper are, for the most part, new 
but it is thought that the generality of the treatment is of value and interest. 


2. General Theory 


The vector V has rectangular components u, v, w which are single-valued and 
differentiable functions of the rectangular co-ordinates x, y, z throughout a certain 
tegion. Two related vectors will be used, namely curl V which has the components 


_ ow 
~ Oy Oz’ Ox’ 


and P with components 


Ou ou Ou 
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Two scalar quantities are employed. These are q, the magnitude of V, given by 
G=u?+v?+w? 


24 oq 
and Q=u ** 


A field line is such that its tangent at any point is parallel to the vector V at the 
point. Then if s is the arc of the field line, we have 


dx: dy:dz:ds=u:v:w:q. ‘ (1) 


It follows from these relations that 


p, (2u , du de) 
“ds ~\ax ds dyds dz ds 
au 
ds * 
Ou ) 
or P, 
dv 
Ow 
Pw 195 


x, Oqdy dq dz 


~ Ox dy ds dz ds’ 
while uP, +0P,+wPe=q(u + wo 


At any point on a field line there is an osculating plane, which is defined as the 
limiting position of the plane through three points on the line when the points tend 
to coincide. That normal to the field line which lies in the osculating plane is called 
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ANALYSIS OF A VECTOR FIELD 


the principal normal while the normal which is perpendicular to the osculating 


plane is called the binormal. We introduce the orthogonal system of direction 
cosines : 


l,,m,,n,, for the tangent 


l,,m,,n., for the principal normal 


for the binormal. 


Let R be the radius of curvature of the field line. 


standard formulae for skew curves : — 


Then we have the following 


1, m, _ d’y n, 
R ds’ R ds?’ R_ ds*° 
Now 


and we obtain by expansion, on account of equations (1), 


ax 


and there are two similar equations obtained by cyclic interchange. Hence we find 


l, =qP,-uQ 
m,(£) =qP,-vQ 
Nn, (4) = qP., - wQ. 


Since the binormal is perpendicular to the tangent and principal normal, we have 


(5). 


1, : Ms :N,=(m.w — N,v): (n,u — 1,w): (1,v 


and, by Lagrange’s identity, 


(m,w n,v)? +(n,u—I,w)? + m.u)? 


=(u? +v? + w’) + m,7+n,7)- (l,u+m,v+nw) 

= 

since |,u+m,v+n,w=0. 
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When we substitute from (5) we get 


3 
(£) =wP, vP, 


m,(£) =uP,-wP, . ‘ (6) 


Let v, be distance measured along the principal normal and let v, refer similarly 
to the binormal. Then 


0q oq 


and, by equations (5), 


+P, 24 | p, <4) - -Q? 


dy 
) oq ) 
=P. (54 -P, +P. -P, 


Multiply equations (5) in order by /,, m,, n, respectively and add: — 


2 


Next, multiply equations (5) in order by P., P,, P» respectively and add: -— 
3 
+ (Pu +P? + Q vPy + 
=q (P.?+P,? + - 
Equation (8) now yields 
a” 
- Q?= ‘ (9) 
Ou dv Ow ou ou ou ) 
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and there are two similar equations obtained by cyclic interchange. By (9) and (10), 
equation (7) becomes, on rearrangement and use of equations (6), 


Ov, 


The expression on the right-hand side of the last equation is the component of 
curl V along the binormal to the field line. 


Next we have —* 


and, by equations (6), 
q oq _ 04 _ =) 4). 
Riv, +P. (wast -uq +Pe (3s 
(12) 
But equations (10) yield 
wa ot =vP,, — wP,+u (ué+vn+ - ‘ (13) 
dy w v. 


and there are two similar equations obtained by cyclic interchange. Accordingly 
equation (12) becomes, on account of (4), 


and finally, by equations (5), 


=~ + man +m. 4) 


The expression on the right-hand side of the last equation is minus the component 
of curl V along the principal normal to the field line. 


When equations (5) are divided by qg they clearly imply that P is the resultant 
of a vector of length qg?/R lying along the principal normal and of a vector of length 
Q=q0q/0s lying along the tangent; this last vector is (0@q/0s)V. Hence we have 
the vector equation 


where N, is a unit vector along the principal normal to the field line. This shows 
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that P lies in the osculating plane. We also deduce from (5) that 


R= P, ux +\ P, +\P, wo): 


In view of equations (2) the last equation becomes 


(16) 


which is an obvious equivalent of a standard formula for the radius of curvature. 


For completeness we add the formula for the tortuosity z of the field line. 


standard formula is 


dx dy dz 
‘ds ds ds 
2 12 
a5 ds 7 as 
d?x 
Now T 73 P. ux 


and it follows that 


Hence by suitable operations on the rows of the determinant we obtain 


u v 
R? P. P, Pw 
OP, OP, dvdq OP» dw dq 


ds ds ds OSs Os Os OS OS OS 


But by equations (2) this becomes 


ueve=w 

R? ow 

tT os -¢s 

au aw 

0s? Os? Os? 


The 


(17) 
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Before leaving the general theory we may remark that when the vector V is 
lamellar, i.e. when it is the gradient of a function of position, its curl is zero, so 
equations (11) and (14) take the very simple forms 


oq 

_ 


3. Applications to Fluid Motion 


Let us first identify V with the instantaneous velocity of the fluid. Then the 
scalar q is the speed of flow, while curl V is the vorticity. In the important case 
where the motion is steady the vector P is now the acceleration of a fluid particle 
at the point considered and we see that the acceleration vector always lies in the 
osculating plane of the streamline at the point considered. The components of the 
acceleration are Q@=40q?/0s along the streamline and q’?/R along the principal 
normal to the streamline, where R is its radius of curvature (and the component 
along the binormal is zero). Now €, 7, ¢ are the components of the vorticity and 
equations (11) and (14) must be interpreted accordingly. It may be noted that these 
are all purely kinematical deductions and therefore true irrespective of the nature 
of the fluid, always provided that it can be treated as a continuum. 


Next suppose that, while the motion is still steady, the fluid is now inviscid and 
barotropic while body forces, if present, possess a potential x. 


Let 6= 
p 
where p is the density of the fluid, and put 
@= - (6+ ). ‘ (20) 
Then the whole set of dynamical equations is equivalent to the vector equation 
P=grad® . : = 


which shows that © is the acceleration potential. If we consider components along 
the streamline, we obtain from the last equation 


a® 
Now let - P=4q?+5+x & 


which is the generalised Bernoulli’s function. 
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Equation (22) can thus be written 


ov 
=0 ‘ ‘ ‘ (24) 


which is the generalised form of Bernoulli’s theorem. 


Now let the fluid be viscous, while still barotropic, and let the motion be steady, 
while the body forces possess a potential. The whole set of dynamical equations 
is equivalent to the vector equation ; 


where v is the kinematic viscosity, and if we take components along the streamline 
we obtain 
dg? 


1 aml; 
(component of VV along streamline). 


This may be rewritten 


~ =v (component of V*V along streamline). . (26) 


If the fluid is incompressible (div V=0) and the flow irrotational (curl V=0) 


we have 


In this case Bernoulli’s theorem (equation (24)) is valid for the viscous fluid. 


As a last application, let us identify the general vector V with the vorticity 
vector, so the field line is now a vortex line. We shall use » for the magnitude of 
the vorticity (in place of q in the general theory) and we shall use the symbol V to 
stand for the velocity. Equation (11) now becomes 


- =(component of curl? V along binormal to vortex line). . (28) 
2 


If the fluid is incompressible we have 


and the equation becomes 
= (component of V?V along binormal to vortex line). (30) 
2 


Equation (14) can be dealt with similarly. 
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Effect of Higher-Harmonic Deflection 
Components on the Creep Buckling 
of Columns 


SHARAD A. PATEL and JOSEPH KEMPNER 


(Polytechnic Institute of Brooklyn) 


SUMMARY: The influence of higher-harmonic deflection components on the 
creep-buckling characteristics of an idealised H-section column is investigated. 
The creep properties of the material of the column are defined by a simple 
power-function creep law. The results show that higher-harmonic deflection 
components reduce the column lifetime significantly only when their initial 
amplitudes, as well as the initial amplitude of the first harmonic component, are 
very large. Furthermore, it is shown that second-harmonic components have 
a much smaller effect on the column behaviour than do third-harmonic 
components. 


1. Introduction 


In recent years several theories of the creep buckling of columns have been 
proposed. In these analyses it is generally assumed that the initial shape of the 
axis of the unloaded column is that of a single half-sine wave and that this shape 
does not change with time (see, for example, Refs. 1-5). In Ref. 5 calculations were 
also made for a non-linear viscous material in which the effects of first and third 
harmonic deflection components were included. 


It is the purpose of the present analysis to investigate the effect of higher- 
harmonic deflection components on the deflection-time characteristics of columns 
whose creep properties can be considered to be non-linearly visco-elastic. 


NOTATION 


A _ total flange area of idealised H-section 
C, = (n?—2)/(1-2) 
E, creep constant 
I moment of inertia of idealised H-section, = Ah?/4 
L_ column length 


gn function of Z,, Z,... 
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distance between flanges of idealised H-section 


h 
m_ creep constant 
n integer 

t time 

ter critical time 

t, | time at which the deflection of the mid-point of the column is h/2 
w deflection due to loads 

w, time-dependent deflection 

w; initial deviation from straightness 


Wr, total deviation of the axis of the loaded column from the x-axis 
at t=0 


x axial co-ordinate 

Zn Fourier coefficient 
Zon = [n?/(n? —2)] Zin 

Zin Fourier coefficient 

€ strain, positive in compression 

creep constant 

o stress 
o, = 7°E,I/(AL’) 

@ average compressive stress 

7 non-dimensional time variable (defined in equation (8)) 
Tx  Critical-time parameter 


7, time parameter corresponding to total mid-point deflection of 
column equal to one-half column depth 


Suffixes 
(1) denotes quantities calculated when z,=0 for n> 1 


x denotes partial differentiation with respect to x 


A dot over a symbol denotes partial differentiation with respect to f. 
A prime (’) denotes final quantities determined by numerical integration. 


2. Differential Equations of Bending of H-Section Columns with 
Initial Deviations from Straightness 
The creep law applied and the basic differential equations derived in the present 
analysis are the same as those of Ref. 4 for a material whose tensile and compressive 
creep properties are identical. 
216 The Aeronautical Quarterly 
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7 
IDEALISED “ONE /2 
7 
=| 
t | a/2 
Fig. 1. Fig. 2. 
Experimental and idealised creep curves. Idealised H-section. 


In the analysis presented in Ref. 4 the creep law was assumed to be of the form 


in which « and co, respectively, are the strain and stress, E,, m, and A are creep 
parameters and the dot over a symbol indicates differentiation with respect to time t. 
In this equation < and o are positive in compression. The creep parameters E,, m 
and A are considered as constants for a given material and temperature. They can 
be obtained experimentally from constant-stress uniaxial creep tests (see Fig. 1). 


The differential equations governing the deflection-time characteristics of an 
idealised H-section column (Fig. 2) which were derived in Ref. 4 are 


(2) 
for W.+ Wro <5 

(3) 


for W.+ Wr => 
in which w, is the time-dependent deflection, accrued for t>0, wr, is the total 
deviation of the axis of the loaded column from the x-axis at time t=0 (Fig. 3), x is 
the axial co-ordinate of a point on the column axis, a subscript x denotes differentia- 
tion with respect to x, A is the cross-sectional area of the column, J/=Ah?/4 is the 
moment inertia of the idealised H-section, h is the distance between flanges, 
and @ is the average compressive stress. Equation (2) is applicable up to the time 
at which the total deflection (w.+w2.) becomes equal to one-half the depth of the 
column; thereafter equation (3) applies. When m is an odd integer, equations (2) 
and (3) are identical. 
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L 
Fig. 3. 
Deflections of 
simply-supported column. 
po 
P 
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3. Deflection of a Column with Initial Deviations 
from Straightness 


Solutions to equations (2) and (3) were obtained in Refs. 4 and 6 with the 
collocation method and the assumption that both the initial deviations from straight- 
ness and the deflections due to loads could be represented by a single half-sine wave. 
However, the column test results reported in Ref. 7 suggest that this assumption 
may not always be accurate, particularly if the amplitudes of higher harmonics are 
of the order of that of the first harmonic. Therefore, in the present work the initial 
deviation from straightness w; is represented by the Fourier series 

NTX 


in which z;, is a constant, L is the length of the column, and the co-ordinate x is 
measured from a support (Fig. 3). Hence, from equation (4) and elastic theory 


n=1 
n? 
CE 
Cy = AL? . . . . (5c) 
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The total deflection can also be assumed as“ 


n=1 


where the Fourier coefficient z, is a function of time only. Elimination of w, and 
We + Wr, from equation (2), with the use of equation (6), yields 


2 


-[1-2 3 asin =0 


When m=1, equation (7) becomes linear and reduces to the corresponding problem 
of a visco-elastic column (see, for example, Ref. (8)). 


The present analysis is confined to a solution corresponding to an odd integral 
value of m, and hence equation (7) is valid even when (w.+ Wr,) => h/2. Because of 
the complexity of equation (7), solutions are obtained only for m=3. This value 
of m was determined from the results of creep-bending tests® on 24S-T4 aluminium 
alloy beams at 600°F. With m=3, equations (7) and (8) become 


sin L sin L +2 Zn SiN L (9) 
# 
tT=2 > t . (10) 


in which C,=(n’?—2)/(1—2). If the right hand side of equation (9) is expanded in 


a Fourier series of the form % g,sin(nzx/L), in which g, is a function of 
n=1 


Zi, Z2, .-.., equation (9) is replaced by an infinite set of first-order non-linear 
differential equations of the form 


dZn 


The corresponding initial conditions are (see equations (5) and (6)) at t=0, Zn=Zyn, 
n=1, 2,... To simplify the calculations it is assumed that only the first three 
harmonics are present initially (z,,=0 for n> 3) and the effects of harmonics 
higher than the third are negligible in the solution (z,=0 and g,=0 form > 3). This 
simplification reduces the problem to the solution of three simultaneous differential 
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equations for the determination of z,, z., and z,. Since a closed form solution of 
these equations could not be found, a numerical method of solution was applied. 
Previous analyses of columns“: » showed that after a sufficient loading time large 
increases in deflection were obtained for small increases in time. Because of this, it 
was deemed expedient to choose z, rather than 7 as the independent variable in the 
present numerical solution. On the basis of the foregoing considerations, the 
resulting equations are 


= (1+ 2z,? + 2,7 + +22,2,) . ‘ 
dz, 1 2 2 2 3 2_37 2 

CD (1+ 2z,? + 2Z,? + 25? —(z, /3) (z,? 32, 
d 


in which D=z, (Z,?— Zo”). 


If z. and z, are taken as zero, equation (14) can be integrated to yield 


1+(1/2)? 


where 7,r;,) is the critical-time parameter corresponding to an indefinitely large value 
of the deflection parameter z, when z,=0 and z,=0. Equations (15) and (16) differ 
slightly from the corresponding relations of Ref. 4 because the collocation method 
was applied in this reference. These results also differ from those of Ref. 5; there 
the elastic part of the creep law (equation (1)) was neglected. Equations (15) and 
(16) are used as a basis for comparison in the determination of the importance of 
the effects of higher harmonics. 


In the present analysis, solutions of equations (12) to (14) are obtained for 
z, 40, z, 40, and z, = 0 and for z, 40, z,0, and z,=0. Thus, specific 
numerical work is done with the assumption that only the first and second harmonics 
are significant and also with the assumption that only the first and third harmonics 
are significant. It may be noted from equation (12) that if z, is initially zero 
(Z.=0) and z,, #0, z, remains zero for all time. On the other hand, equation (13) 
reveals that if z,,=0 and z,, 0, z, will develop with time. These observations 
are consistent with the physical problem. Thus, because of the symmetry of load 
and column about the point x=L/2, odd-harmonic initial deviations can result in 
only (but all) odd-harmonic deflections. 
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When z, =0, equations (13) to (14) yield 


dz, _ ((1—2)/(4—2)] z, (1+2z,?+z,”) 


dz, Z, (I +2,? +22z,”) 
dr 2 
dz, 3z, (1+z,?.+2z,) (18) 
with the initial conditions z,=z,, and z,=Z). at r=0. 
When z, = 0, 
dz, Z, (1+z,? + 2z,?— 2,25) 
dr y (20) 


dz, 3z, 


with the initial conditions z,=z,, and z,=Z,, at *=0. 

The two sets of equations were solved by the Runge-Kutta method". Thus, 
for example, z, was determined numerically from equation (17) as a function of z,; 
thereafter + was determined as a function of z, from equation (18). In the latter 
work the Runge-Kutta method reduces to Simpson’s rule. Similarly, z, and + were 
determined as functions of z, from equations (19) and (20). Calculations were 
performed for the following values of the parameters : — 


Zo1 =0:01; (Zi: =0-003) 


1 


1 
Zis/Zu=- -1, 
and 
Zo, =0°1; (zi, =0-03) 
Zio / Zi ‘3 
a=0-7, i2 = 4 

1 


The corresponding values of z,, and z,,, computed with the aid of equation (5a), are 
given in Table I. The values of the ratio z;,/z,, were chosen in accordance with the 
results of measurements reported in Table I of Ref. 7. It may be noted from 
equations (17) and (18) that, if a set of values z,, z,, 7 is a solution, then the set 
Z,, — 2, 7 is also a solution. Thus, it is only necessary to investigate positive values 
of the ratio z;,/z;,. However, a similar condition does not exist for equations (19) 
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Fig. 4. 


Variation of amplitudes of first and third harmonic deflection components. 


and (20). The calculations were confined to negative values of the ratio z;,/2i,, since 
these correspond to larger initial deviations of the mid-point of the column than do 
the positive values and, hence, to smaller values of the critical time. 


In each case considered, the computations showed that z, increased much more 
rapidly with time than either z, or z,. Furthermore, d7/dz, appeared to approach 
zero with increasing z, in the manner indicated in Fig. 4. Hence, the calculations 
were performed as described previously until a sufficient time 7’ was obtained such 
that the corresponding displacement amplitude z,’ was very much greater than either 
Z.=Z,’ OF Z;=2Z,'.. Thus, equations (18) or (20) reduce to 


2,’ 
or r=1'+ alo 5 122) 


From equation (22) it is seen that as z, grows indefinitely large s approaches the limit 
, 
Top = + 3 log (1+(1/z,)?] . 


In all the calculations z,’ is taken equal to 3, and therefore equation (23) becomes 


Examination of equations (17) or (19) for z, >> z. or z, >> Z, in the region 
Z, =z,’ shows that z, and z, remain much smaller than z,, although they increase 
indefinitely. with z,. 
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TABLE 1 
EFFECT OF SECOND- AND THIRD-HARMONIC DEFLECTION COMPONENTS ON CREEP BUCKLING 


(a) z, 40, z,0, z,=0; a=0°7 
Zo, = 9°01; z,=0-003; z,’=3 


Tera) = 3°07; =2'53 
090909 16955 3-037 3-07 2:53 1 1 
1 3°6364 6°7827 3-037 3:07 2°53 1 1 
3 10-909 20°346 3:037 3-07 | 1 
Zy,=0'1; %,=0-03; z,’=3 
= 1545 74 (y= 100 
4 9-0909 13-746 1-504 1:54 1-00 1 1 
1 36°364 53°811 1:500 1°54 1:00 1 1 
3 109-09 164-20 1°443 1-48 0:952 0:96 0:95 
(b) z,=0; a=0°7 
Zp, =O01; =0°003; z,’=3 
= 3-075 744) =2'53 
j —~Zo3 x 10° x 10% r Tor For 
081325 3-036 3-07 2:53 1 1 
1 3-2530 25-853 3-034 3-07 2:53 1 1 
3 9-7590 34-129 3-030 3-07 2:52 1 1 
) = Ol; %,=0°03; z,’=3 
= 1545 = 100 
) 0 0 2:2685 1-503 1°54 1-00 1 1 
4 8-1325 31:349 1502 1-54 0999 1 0:99 
1 32-530 59°638 1-485 1:52 0-950 0:99 0-95 
3 97:590 13484 147 0841 095 0:84 
3) For each case considered, Table I shows the corresponding values of z,’, 7’, and 
ty OF Z,’, 7 and t.,. Also included are the values of the time parameter 7, at which L 
es the mid-point of the column exhibits a total deflection equal to one-half the column 
depth. For equations (17) and (18), this corresponds to z,=4; for equations (19) 
4) and (20), to z, -Z,=4. When only first-harmonic solutions are considered, equation 
(15) applies with z,==4 and yields 
on 
se 


erly August 1957 223 


| 


SHARAD A. PATEL AND JOSEPH KEMPNER 


The parameter 7, corresponds to the time at which the stress in the flange on the 
convex surface of the beam vanishes; thereafter it becomes tensile. Earlier calcula- 
tions indicated that the deflections increase very rapidly after 7=7, (Ref. 4). Also, 
in Ref. 9 it is suggested that column lifetime be characterised by a finite limiting 
deflection sufficiently small so that the small-deflection bending theory used in the ( 
analysis is justified. Reasonable correlation between theory and experiment was f 
achieved when the collapse time of the column was taken to correspond to 7,. 


4. Results and Discussion 


The last two columns of Table I present the ratio of the critical time f,, 
obtained when higher harmonics are considered to the critical time t..;,) corres- 
ponding to the consideration of only the first harmonic equation (16), as well as the 
ratio t,/f,¢.)=71/71,). These results show that, when the amplitude of the first 9 
harmonic of the initial deviation is small (z;,=0-003, z,,=0-01), the presence of 
higher harmonics, even of the same magnitude as the first harmonic, has negligible 
effect on the lifetime of the column. On the other hand, if the initial deviation is 
large (z;,=0-03, z,,=0°1), the effect of higher harmonics may be appreciable, 
particularly when the amplitude z,, of the third harmonic is of the same order of 
magnitude as the amplitude z,,. In the latter case ¢t,/t,,,,=0°84. The results 
also indicate that third-harmonic deflection: components have more effect on the 
column life than do second-harmonic components. 


Since only unusually large initial amplitudes of higher-harmonic deflection 
components influence the column lifetime significantly, it is sufficiently accurate to 
consider only first-harmonic amplitudes in analyses of creep buckling of columns. 
It is, however, important that an accurate estimate of the first harmonic amplitude, 
rather than simply the total deviation of the mid-point of the column be used. 
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The Oscillating Two-Dimensional Aerofoil 
Between Porous Walls 


D. G. DRAKE 


(Department of Aeronautical Engineering, University of Southampton) 


SUMMARY: The compressible flow past an oscillating two-dimensional aerofoil 
in a wind tunnel with porous walls is considered, using linearised theory. The 
porous wall is assumed to have the property that the ratio of the normal velocity 
at the wall to the pressure drop across the wall is constant. Transform theory is 
used to find the supersonic longitudinal stability derivatives, and an extension 
of Possio’s integral equation for the quasi-stationary case in subsonic flow. 


1. Introduction 


The possible advantages of a wind tunnel with porous walls, compared with a 
solid-wall tunnel, have been discussed by Baldwin, Turner and Knechtel™, Woods”? 
and others. Although the blockage factor may be made to vanish when the walls 
of a wind tunnel are porous, wall interference effects still exist and must be evalua- 
ted in order to correct wind tunnel measurements to free air conditions. Expressions 
for the supersonic longitudinal stability derivatives are found in Section 4. 


One method of representing the effect of solid tunnel walls and free jet 
boundaries is by an infinite set of images, but it is difficult to see how to apply this 
method for the case of porous walls. The transform technique, employed by 
Miles® for the solid-wall case and by the present author“) for the free-jet case, 
enables the boundary condition at the porous walls to be satisfied. 


Consider a uniform stream of velocity U flowing past an aerofoil of chord c 
which lies midway between two porous walls at a distance b apart. Let the aero- 
foil perform simple harmonic oscillations of complex amplitude — U2 (x) with an 
angular frequency » such that the normal velocity component is -U & [a (x) e’'], 
where & denotes “ real part of”. The co-ordinate system, moving with a velocity 
U relative to the main stream, consists of x- and y-axes, parallel and normal to the 
walls respectively, with origin at the leading edge of the aerofoil. The restriction 
| «| << 1 allows the linearised theory to be used and the pressure p, after making 
x and y co-ordinates non-dimensional with reference to c, is given by 


p(x. +ik) (x, y) : (1) 
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where &["/coe] is the perturbation velocity potential, k=wc/U is the reduced 
frequency and p, is the density of the undisturbed flow. Since ¢ is anti- 
symmetrical in y, the pressure jump across the aerofoil is given by (1) as 


where y(x)= (> +ik) (x, 0+) (3) 


If n is the co-ordinate along the outward normal at the walls, then the normal 
velocity at the walls is % [U¢,e]. The porous wall is assumed to exhibit the 
property that the normal velocity at the wall, #[U¢,e], is proportional to the 
pressure drop Ap across it, the constant of proportionality being p,U/o, with o 
known as the “ porosity ”, i.e. 


2 
One of the equations of motion is 


PU? $2) e]=p,—p, « « 


where p, is the free stream pressure. With the further assumption that the pressure 
outside the wall is equal to the free stream pressure, the equation of motion (5) and 
equation (4) lead to the boundary condition 


Therefore 


%)=0, y='+b/(2c) =ity,, say. . 


The boundary condition (4) is used by Baldwin, Turner and Knechtel®, Woods’? 
and Goodman™. The special cases of the porosity being zero and tending to 
infinity give the solid-wall tunnel and the free jet respectively. 


NOTATION 
x,y non-dimensional Cartesian co-ordinates 
chord of the aerofoil 
b distance between the porous walls 


free stream density 
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U free stream velocity 
» angular frequency of oscillation of the aerofoil 
k = wc/U, reduced frequency parameter 
n co-ordinate along outward normal at the walls 
M free stream Mach number 
p pressure 
Pp, free stream pressure 
t time 
KR “real part of” 
@ velocity potential 
Laplace transform of 
a complex amplitude of incidence 
a Laplace transform of 
y defined by equation (3) 
y Laplace transform of y 
g defined by equation (17) 
g inverse Laplace transform of g 
porosity 
o’ = 2cy,, porosity parameter 
y, = b/(2c), half the non-dimensional distance between the porous walls 
B? = 
s Laplace transform variable 
f(s) Laplace transform of f (x) 
A = B[(s+iKM) + K*} 
K = kM/B?= —kM/f? 
vy Fourier transform variable 
F(v) Fourier transform of f (x) 
A = B[-(v+KM)?+K?}! 
A = 2By,, tunnel aspect ratio 
I’ Fourier transform of y 


® Fourier transform of ¢ 
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defined by equation (27) 
inverse Fourier transform of ®, 
“inverse Laplace transform of” 


lift and moment coefficients respectively 

defined by equations (37) and (38) 

co-ordinate of the centre of pressure 

point about which pitching motion occurs 
amplitude of pitching motion 

positive integer determined by N<1/A<N+1 


1 (B-c) 
tos. (Bto) +2nei | 


yi 


denote partial differentiation 


2. Supersonic Free Stream 


For irrotational motion the equation of continuity is 


— $,,=0 


so that 


py (x), y=0. 


The boundary conditions at the porous walls are 


oy y=+y 


and the supersonic flow initial conditions 


¢=0.x=0. 


and ¢,=¢,=0, x<0. 


Using the Laplace transform defined by 


f(s)= | e-f (x) dx, 


equation (8), together with (11) and (12), gives 


229 


(8) 


The use of the linearised theory allows the downwash condition at the aerofoil to be 


0) 


(10) 


(11) 
(12) 
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where \?= [(s+iKM)?+K?] and K=kM/B?’. 
Under this transform, (9) and (10) become 

and (s+ik) ¢=0, y=H+y,. . 


The solution of (13) subject to (14) and (15) is 


cosh A (y, —| y |)+o (s+ik) sinh A (y, y sgn y 


[> sinh (Ay,)+o (s+ik) cosh ay.) | 


[a cosh A (y, — | y |) +o (s+ik) sinh A (y,- | y sen y 
where g2(s, y)= 


(17) 
sinh (8+ ik) cosh ay.) | 


An expression for y(x) is needed when determining the lift and moment 
coefficients of the aerofoil. Substituting from (16) into (3) and utilising the Faltung 
product theorem 


x 
(x)= +ik) | g(x-z,0+)a(z)dz. . . (18) 
0 
From a theorem in transform theory 


Lim g (x,0+)=Lim sg (s,0+)=B-'. (19) 
x0 


and, therefore, by differentiation 


y (x)= (x)+ h(x—2z) «(z) dz, 


where (= + ik) g(x, 0+). 


3. Subsonic Free Stream 


In subsonic flow the conditions downstream of the aerofoil influence the condi- 
tions at the aerofoil. The downwash is unknown away from the aerofoil along 
y=0, since the potential is discontinuous across the wake, and an alternative 
condition on the potential must be sought. 
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A convenient condition is on the pressure, which must be continuous except 
across the aerofoil. Since ¢ is anti-symetrical in y, equation (1) gives 
¢,+iko=0, y=0,x<0 and x>1. (20) 


The formulation of the subsonic problem is, therefore, that ¢ satisfies 


with the conditions dyto (o2+¢)=0, y=+y, . (22) 
and drt +ty(x), y 20. (23) 


Using the Fourier transform F (v)= | (x) dx,  (21)-(23) become 


+o (iv+ik)®=0, y=+Hy, ‘ (25) 
and . . (26) 


where A?=B?[-(v+KM)?+K?] with K=- kM/B’, 
and I" (¥) is the Fourier transform of y (x). 


The solution of (24) subject to (25) and (26) is 


cosh A (y, | y |)+io (v +k) sinh A (y,- | y | )| sgn y 


y)= 
i(v+k) [a cosh (Ay,)+io (v+k) sinh (Ay,) | 


[ A cosh A (y,—| y D+ io +4) sind (y, y sen y 
where ®, (v, y)= 


i(v+k) [:a cosh (Ay,)+ic (v +k) sinh (Ay, | 


Applying the product theorem (27) gives 


1 


@ 
-Lim 6, @-z @ dz, (28) 


where 9, is the inverse Fourier transform of ®, and the range of integration is 
restricted by (20). 


August 1957 231 


| 
) 
| 
) 
t 
) 


D G. DRAKE 


The inverse of ®, may be found by evaluating the integral 


co 
¢, (x, y= D, (v, y) dy, 


but for general k the work is complicated. If k is assumed small the integral may 
be evaluated by an argument similar to that of Miles“. The assumption is made 
that k has a small negative imaginary part to satisfy the null condition at infinity, 
and later letting the imaginary part of k tend to zero. With this assumption the 
contour of integration in the v-plane may be taken as the real axis with a semi-circle 
of infinite radius in the upper or lower half planes according as x >0 or x <0. 
Evaluating the residues at the poles of ®, exp (ivx) and neglecting 0 (k*) terms, but 
treating ok as small only when compared with o leads to 


[8?-cik (y, - Ly DI 
(x, y)= e sgn y ikoy,}? 
e 


(1 + e~ ‘KMx sin (b,y) 


bay, b,B 


Ly 005 (buy) sin ay) (4 + 5) ] 


where b,y,=tan-! (8/c)+nr, with and x>0. 


When x<0, 
(x, y= (1 + sin (a,y).+ 
n=0 Gn; 
ioke-nx/8 2 
(62 +0") [» cos (any) — sin (apy) (3 + =)] (29) 
where with <0. 


Higher order terms may be included in the k-series expansion of ®, but the 
work involved in the inversion of such terms becomes excessive. Woods tackled 
the steady problem and a solution of the unsteady problem including terms 0 (k) 
can be found from it by iteration. 


Equations (28) and (29) are an extension of Possio’s integral equation. It may 
be seen that the limiting cases of s —> 0 and « —> © correspond to the solid-wall 
tunnel and free-jet™) quasi-stationary results. 


232 


The Aeronautical Quarterly 


ra 


he 
| 


OSCILLATING AEROFOIL 


4. Supersonic Lift and Moment Coefficients 


The supersonic lift and moment coefficients may be found, to the first order of 
k, for a particular prescribed motion of the aerofoil. For a pitching motion of small 
amplitude «, about x=x,, 


a (x)=, [1 +ik (x- x,)] 
and «@(s)=«, [s-!+ik (s-?-x,s~')]. 
Also, the transform of (3) gives 
The lift coefficient is given by 


x 


Since the Laplace transform of { f (x) dx is sof (s), 


where ~-* denotes “inverse Laplace transform of.” Similarly, the moment 
coefficient 


(x, x) y (x) dx 


0 
Neglecting O (k?) terms, A= B (s+iKM) and, by the shift theorem, 


Bs cosh (Bsy,).+ (s— iKM +ik) sinh (Bsy,) 
Bs [Bs sinh (Bsy,)+o (s—iKM + ik) cosh (Bsy,)] 


g(x,0+)=(1-iKMx) { 


=(1-iKMx) { B cosh (Bsy,)+o sinh (Bsy,) 


Bs [B sinh (Bsy,) +o cosh (Bsy, 


= oik 
+2 { B’s? [B sinh (Bsy,) +o cosh (Bsy,)} } 


to the first order of k. 
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c+ico 


The inverse of f (s) is given by 4 { e“f (s) ds, 


c-ioo 


where c is a real number, and the inverse of the expressions in (34) may be found by 
evaluating this integral. The contour is taken to be an infinite straight line parallel 
to the imaginary axis in the s-plane, a distance c from the origin such that c is 
greater than the real part of the positions of any of the poles. The contour is 
completed by an arc of the circle with centre the origin and radius tending to infinity 
on the negative real axis side of the infinite line. 


The zeros of B sinh (Bsy,)+o cosh (Bsy,) are at 


(B- 
(B+ 


where A=2By, is the “tunnel aspect ratio”. The significance of A is that when 
A > 1, i.e. (M? - 1)! b >, the flow at the aerofoil is not influenced by the presence 
of the porous walls. 


Then (34) becomes, with s, given by (35), 


1 iKMx  2iky, ikx 


45 
n= —co By,s, * 5,7Bty,? (B?-o?) s,°Bty,? (B? B*sny, 


Using the Faltung product theorem, equations (30) - (32) and (36) give 


Mt 2ik 
Ci =42, 3 { (1 —ikx,)+2ik } x - > 
ikt | cikte’" } dt, 
+ By s, * (Bo) B’s,y, 


evaluated at x=1. Integrating the infinite series term by term 


C,=42, ikx, +ik— aM + 


1 1 
by (ikx - 1 - 2ik)- 


B*y Yi n=-co 


oik 3cik co 
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co 
The infinite series >» 


1 _ 
4c? 
eo 3 2_ R2 
g 1 _ 


Sn 


1 At (co? At (c? B?) 
oo n 
The infinite series <>. 


are Fourier series and it may be shown that 


1 B B?) (co? B’y = B?) 


where, in each case, N is a positive integer given by N<.1/A<N+l. 


A partial check may be applied to these results by considering the limiting cases 
ao —> 0 and « —> 00 when the expressions for the solid-wall tunnel® (¢ —> 0) and 
the free jet“ (¢ —> 00) are recovered. The expression for C; at «=0 is indeter- 
minate, which is due to the order of the pole at s=0 in the terms of (34) changing 
when o=0. However, the limit as o—>0O does exist and gives the solid-wall 
results. This limit may be found by repeated application of L’Hospital’s rule to the 
coefficients of the independent variables in C;. The limit as « —> 00 is straight- 
forward. 
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Fig. 1. 
The curves of (B/4) I against A for various values of o’. 
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Position of the centre of pressure x, plotted against A for different values of o’. 
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Expressing the lift coefficient in the form 
Ci =, +0 (k*)], . . (37) 


it may be seen that (B/4) J. is a function of the “ tunnel aspect ratio” A=2By, and 
the “ porosity parameter” o’=2cy, only. In Fig. 1 a set of curves is obtained by 
plotting (B/4)/. against A for various values of the parameter o’. The curves of 
(B/4) 1. for non-zero porosity are similar to the free-jet curve and differ from the 
solid-wall curve in so far as (B/4)l1.—>0 as A-—>0O for o’ +0, whereas 
(B/4)1.—> 00 as A —> 0 for the solid wall. The curves of (B/4) 1. have discontinui- 
ties in slope when 1/A=N, N=1,2,...... , and for the special cases of the solid- 
wall tunnel and free jet the slopes are constant between these discontinuities. As 
already stated, when A > 1 the porous walls do not affect the flow at the aerofoil 
and the value of /. is 4/B, independent of the porosity parameter o”’. 


It is seen that, knowing the porosity of the tunnel walls, the tunnel dimensions 
and the Mach number, it is possible from such a curve to adjust an experimentally 
determined /, to free-air conditions. The dependence of /, on the porosity, tunnel 
height and Mach number may be exhibited in a similar manner. 


The moment coefficient may be found from (30), (31), (33) and (36) in a similar 
manner, and may be written in the form 


Cu= 44, {+ (2ikx, $y — + 


+ _ —ikey—1) + (ik+ikM-%) + 
30 


52 2ikx, —ikx,? -ik+ x, By. 

B’y, Sa" B’y, (B? - Sn 

2cik 


where the sums of the infinite series are as in the lift coefficient. Writing 
Cn=%, [m. + ikma+ O (k?)], ‘ (38) 


the position of the centre of pressure is determined by equating m. to zero. 
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This gives 

C.D. 

a. Fi, 

a 

o n=-a Sn n=-co Sn 


Figure 2 shows the position of the centre of pressure x.» plotted against A, with o’ 
as the parameter. 


These curves may be used in the same way as those of Fig. 1 to estimate porous 
wind tunnel wall effects on experimentally determined centres of pressure. For 
non-zero porosity it may be shown that in the limit as A —> 0 the position of the 
centre of pressure x.,, Moves in an oscillatory fashion between the limiting values 
x=4 and x=(l-o’)/2. In the special case of the solid-wall tunnel 


For the free jet the centre of pressure oscillates between x=4 and x= -0O and 
never moves aft of the mid-chord point. 
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An Experimental Investigation of the 
Secondary Flow Occurring in a Compressor 
Cascade 


W. D. ARMSTRONG 
(Cambridge University) 


SUMMARY: When a non-uniform stream flows through a cascade of compressor 
blades a three-dimensional secondary motion is generated in the outlet stream. 
This motion, which consists of a spanwise re-orientation of the Bernoulli planes, 
together with the effects of vorticity components in the stream direction, may 
be predicted theoretically. In this paper experimental results obtained on a 
6 in. chord cascade are reported. The experimental! and computed outlet angles 
are compared and the discrepancies are explained qualitatively. Additional 
particle track and visualisation experiments were designed to investigate the 
local blade flows and the effects of the nose vortices. 


1. Introduction 


In recent papers“: ” it has been shown that ideally the secondary flow which 
occurs when a non-uniform stream passes through a cascade of turning vanes or 
blades can be computed by following the movement of the vorticity vectors. When- 
ever the upstream vorticity or the angle of turn is small, the perturbation method of 
Squire and Winter or Hawthorne“ may be used for the calculation of the down- 
stream vorticity component in the stream direction between the blades. 
Qualitatively this vorticity component may be imagined in terms of induced velocity 
components. The pattern of these for a wall boundary layer flowing through a 
cascade of blades is shown in Fig. 1. 


Experiments have been devised to measure the induced velocities and thereby 
compute the secondary flow, but it was found at an early stage that in many cases 
of practical importance a simple model such as shown in Fig. 1 could not be used to 
visualise the flow. Another paper“, however, gives details of tests using an impulse 
turbine cascade in which the simple theory could be applied. 


In the experiments presently reported, in addition to velocity measurements the 
three-dimensional tracks of air particles have been constructed, so enabling a model 
of the flow pattern to be built. 
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SECONDARY FLOW IN A COMPRESSOR. CASCADE 


VELOCITY 
DECREASING 


Fig. 1. 
Induced velocity components. 


NOTATION 

a, air inlet angle 

@ air outlet angle 

Q inlet vorticity 

secondary stream function 

¢ stagger angle 

X, 
X.Y.Z | rectangulas co-ordinates (see Fig. 2) 


U, upstream velocity 


2. Experimental Details 
2.1. THE APPARATUS USED 


A cascade of nine compressor blades was used in a 150 h.p. pressure wind 
tunnel. The blade profile, C4, was set on a 30° circular arc camber line with a 
chord of 6 in. The blade spacing was also 6 in. and the aspect ratio 3. A stagger 
angle of - 36° was chosen in conjunction with an air inlet angle «, of 50° to give 
a design air outlet angle of 30°. This blading is shown in Fig. 2, which also 
indicates the position of a boundary layer suction slot ahead of the convex end 
blade of the cascade. The suction was adjusted to give a uniform flow upstream of 
the cascade along almost its whole length. Downstream of the blade trailing edges 
the walls of the cascade were extended for a distance of 14 chords. 


Normally the boundary layer on the side wall of the tunnel which passes over 
the blade ends has a thickness of 1 in. at the Reynolds number used (3-9 x 10°), but 
for these tests it was artificially thickened to give the spanwise upstream velocity 
profile shown in Fig. 3. The artificial layer was created by increasing the length of 
the parallel portion of the tunnel and inserting within it two sections of bakelised 
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3FT.3IN.TO VELOCITY 
PROFILE HONEYCOMB 


END OF - 
OOWNSTREAM 
WALLS 


AIR DIRECTION 


SCALE }—SIN-_ 
CHORD GIN. PITCH GIN. SPAN IBIN.o< = 50°C=-36° 
1oca/socso 


40 
Fig. 2. 


Diagrammatic arrangement of compressor cascade. 


paper honeycomb of triangular plan form—the base of the triangle being along the 
wall and the apex pointing towards mid-stream. This formed the foundation for 
strips of stainless steel wire gauze which were added on the upstream side of the 
honeycomb to increase its effect. 


2.2. AIR FLOW MEASUREMENTS 


The instrumentation developed for this work has been described elsewhere“; it 
enabled traverses to be made one chord upstream of the leading edges and in any 
plane downstream of the trailing edges. Observations, near the trailing edge and 
one chord downstream, were made of the total head and the two angles, which 
enabled all the velocity components to be computed. 


2.3. PARTICLE TRACKS. TESTS WITH HoT WIRE AND THERMOCOUPLE 


As a convenient means of locating the tracks downstream of the blade trailing 
edges, a small filament of air was continuously heated in its passage through the 
blading and then the warmed air was located in a parallel plane farther downstream 
The filament of air was heated by a short length, approximately ? in., of 28 s.w.g. 
resistance wire which was maintained with a red glow by a current of up to 10 amp. 
The resistance wire was suspended on two legs of copper wire through which the 
heating current was supplied. Several leg lengths were used (between } in. and 
$ in.) and the bases of the legs were fixed to, various positions on the blade trailing 
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SECONDARY FLOW IN A COMPRESSOR CASCADE 


0-9}- 
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U 
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° 
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J 


| CENTRE LINE 


— L L 
WALL 2IN. 4iN 6IN. SIN. 
Fig. 3. 


Upstream velocity profiles in the centre of the blade passage with artificially thickened 
boundary layer. 


edge or the cascade wall. The copper-constantan thermocouple used to locate the. 
heated filament of air was traversed in the normal yawmeter traversing gear in a 
plane one chord downstream of the trailing edges. A balancing thermocouple was 
located in the same plane, in the air flow from another blade passage, and the 
difference between the temperature recorded by the pair of thermocouples was 
measured by a swinging coil galvanometer with large optical magnification. For 
each position of the heating element the traversing thermocouple located a small 
region of heated air but, by cross plotting, the position of the peak air temperature 
could be accurately located. 


2.4. Low SPEED VISUALISATION EXPERIMENTS 


For the detail investigation of the flow within the boundary layer a low speed 
(5 ft./sec.) open tunnel discharging over a flat plate was used. Smoke filaments 
were injected into the boundary layer a short distance along the plate, through a 
passage whose open end was flush with the surface. Models of blades were stood 
on the plate in the boundary layer stream and photographs of the smoke filaments 
were taken. 
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Fig. 4. 
Total head in the plane of the trailing edge expressed as a fraction of the mean centre-line 
value. 


3. Results 


3.1. THe ToraL HEAD DISTRIBUTION DOWNSTREAM 


Complete traverses covering 1000 points were made over half the blade span 
and two passages, in each of the two downstream planes. The total head contours 
for these two positions are plotted in Figs. 4 and 5. On these maps the conventional 
secondary flow—a pair of vortices rotating near each trailing edge and a pair of 
vortices within the blade passages having a rotational direction opposite to that of 
the trailing wake vortices—is clearly seen. In addition there is a large region in the 
corner of the blade passage, near the blade suction surface, where the total head is 
zero. Within these regions the flow appeared to be very unstable and it was often 
impossible to balance the yawmeter, the flow apparently being reversed. 
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Total head one chord downstream of the trailing edge expressed as a fraction of the mean 
centre-line value. 


3.2. THE OUTLET ANGLE 2, 


In the plane of the trailing edge the outlet angle varies widely, both across a 
pitch and across the span. Because the region of unsteady flow in the corner 
between the suction surface and the wall prevented accurate measurement of 2, 
there, it is not possible to present integrated values of ~, along the span. The 
variation along a line midway between the trailing edges is plotted in Fig. 6. Fora 
position one chord downstream, two representative curves are shown in Fig. 7. One 
curve gives the variation of <, along a line parallel to the blade trailing edges and 
midway between the blades, while the second relates to positions along a curved 
line lying midway between the wakes, the positions of which were determined in 
the experiment. 
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Fig. 8. 
Secondary flow stream function v. 


(Ince UNITS) 


Measurements of the inlet vorticity enable the secondary vorticity to be 
calculated by using the theory of Hawthorne“ and assuming that the angle changes 
linearly with the distance along the chord. The distribution of the secondary stream 
function y in a plane perpendicular to the main stream direction may then be com- 
puted by relaxation techniques. This theoretical distribution of ¥ is shown in 
Fig. 8. From this it will be seen that the induced velocity perpendicular to the 
span, 0 /0z, should become zero at about 3-3 in. from the cascade wall, i.e., in the 
plane of the trailing edge, z, at 3-3 in. from the wall and in the middle of the stream 
should be the same as the two-dimensional centre-line value, to a first approxima- 
tion. The distribution of the theoretical z, is shown on Fig. 6. There are obvious 
discrepancies. 

August 1957 247 


BLADE 
j. CENTRE 
TIE. LINE 
| BIN. 
0-10 | 
| | 7IN. 
| 0-20 | 
6IN. 
| 0.30 | 
| 
0-40 | 
| 
| | 
2 IN. 
| | LIN, | 
= 
| 
| 


Ww. D. ARMSTRONG 


| 
t 
' 
\ 
\ 
\ 
4 
38 
“, 
36] 
=) 
EXPERIMENTAL 
z 
° > 
32 
» 
ay 
WALL 2iN, 4iN. 6IN. BIN. 
Fig. 9. 


Outlet angle «, one chord downstream of the trailing edge at the mid-wake position. 


Between a plane including the trailing edges and a plane one chord downstream 
there are very considerable changes in «, (compare Figs. 6 and 7). These 
re-orientations cannot be due to the secondary rotational flow, but are apparently 
the result of two interlinked effects, local stalling of a portion of the blade and 
readjustment of the axial velocity profile according to the actuator plane theory™. 
The former effect is discussed in connection with the particle track experiments in 
Section 3.4. 


3.3. THE APPLICATION OF THE ACTUATOR PLANE THEORY 


It has been shown™ that, even in the absence of secondary vorticity components 
in the stream direction, there is a change in the air outlet angle downstream of a 
blade row when the upstream flow is non-uniform across the blade span. This is 
necessitated by the equations of motion and continuity, which can be linearised and 
solved if the blade row is replaced by an “actuator plane”. To demonstrate the 
direction of the a, variation, the theoretical values of «, at an infinite distance 
downstream have been calculated using the actuator plane theory; they are shown 
in Fig. 9, along with the experimental inter-wake values one chord downstream. 
(a, at the actuator plane was assumed constant). It will be seen that the actuator 
plane theory predicts changes in «, of the same order as those measured. 
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Fig. 11. 
Photographs of goosedown attached to the cascade wall and blade. 
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Fig. 12. 
Cascade flow in the very low speed tunnel. 
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Fig. 13. 
The nose vortex on a large single aerofoil. 
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3.4. THE PARTICLE TRACK EXPERIMENTS 


The results are presented in the form of a model, Figs. 10(a) (6) (©). In 
Fig. 10(5), which looks down onto the concave surfaces of the blades, the left hand 
side represents the cascade wall and the right hand side the centre line of the blade 
span. The blades themselves may be seen in the other photographs; Fig. 10(a) is a 
view through the cascade wall towards the blade centre line, while Fig. 10(c) is 
taken from the centre plane. 


In the blade passage the model clearly shows the secondary rotational flow; the 
air leaving the trailing edges has observable components of velocity in the direction 
of the blade length. These are towards the blade centre line on the convex surface 
and towards the wall on the concave surface. 


Within the wall boundary layer region, the particle track experiments confirm 
conclusively the results obtained in the total head and angle traverses. Near the 
blade suction surface and the passage wall was found a region where the heated air 
was dispersed throughout the flow. This region corresponds to the zero total head 
area. The change in «, in passing from the trailing edge to a position one chord 
downstream is clearly shown in the photograph of the model looking through the 
side wall (Fig. 10(a)). It is apparent that any continuation of the overturning 
(i.e. reduced a) which occurs within the blades due to the secondary flow is opposed. 
by the tangential movement of the high energy air which has come along the concave 
surface. This higher energy air moves into the region near the convex blade 
surface where the total pressure is low and then spreads out along the cascade of 
blades. It gives the effect of a continuous underturning of the air all along the blade 
span. In this case, after a distance of one chord there is little or no air deflection in 
the air flowing near the cascade wall, although within the confines of the blades the 
overturning amounts to several degrees. 


It is interesting but experimentally difficult to trace the movement of the region 
of low total head—the region near the trailing edge which appears to have collected 
most of the low energy flow from the whole of ihe cascade side wall. Apparently 
this air slowly becomes dispersed some distance away from the wall in a pocket 
which extends the wake width. If boundary layer suction on the annular walls of a 
compressor were a practical proposition, it appears that the most convenient place 
to apply the suction would be in the corner between the blade suction surface and 
the wall, where the secondary rotational flow has collected together all the low 
energy air. 


Photographs of goosedown attached to the cascade wall and blade are included 
in Fig. 11. These also show the stalling of the downstream part of the blade suction 
surface or the separation of the blade boundary layer from this surface due to 
transport of the wall boundary layer onto the blade and to the change in outlet 
angle predicted by the actuator plane theory. 


Experiments with these compressor blades have been carried out at other 
incidences and with thinner cascade wall boundary layers, but in all cases a region 
of stalled flow was observed in the corner. The extent of the disturbance, however, 
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is greater with the thicker boundary layers and at the higher incidences. This effect 
may also be the primary cause of the high contraction coefficients observed in 
cascade tests at high incidences. 


3.5. MODEL EXPERIMENTS IN THE VERY LOW SPEED TUNNEL 


A cascade of 3 in. chord compressor blades of the same profile but with a 
pitch/chord ratio of 1/2, was also set up in the low speed smoke tunnel. Two 
typical photographs of the flow observed in the cascade are reproduced in Fig. 12(a) 
and (5). In (5) there are two smoke traces; one lies within the boundary layer on 
the plate and is emitted from the hole near the leading edge of the plate, while the 
other emanates from the trailing edge of a laminar flow aerofoil held in the main 
stream. This latter smoke filament indicates therefore the turning of the main air 
stream. The simultaneous smoke traces serve to show the effects of secondary flow. 
Just as in the large cascade, the flow on the blade suction surfaces near the trailing 
edge has separated and the smoke filament usually appeared to “corkscrew” 
upwards out of the air stream and along the trailing edge of the blade. Downstream 
of the trailing edge it is just possible to distinguish the change in air deflection of 
the main stream coming over the pressure surface. 


Model tests in this tunnel, using vortex generators upstream of the leading 
edges, showed that it was possible to modify the upstream velocity profile 
sufficiently to reduce the main component of distributed secondary vorticity. This 
reduced the apparent stalling, but full scale trials using the 6 in. chord cascade 
were not so successful. 


In the smoke filament experiments some interesting effects were observed near 
the rounded leading edge of each aerofoil of the cascade. This flow was investiga- 
ted in greater detail by using a larger single aerofoil. The photographs of Fig. 13 
show what is now called the “nose vortex”. This is a large well-defined rolled-up 
vortex or group of vortices which passes round the nose and trails off with the main 
stream on each side of the aerofoil. In addition to the single one shown, other 
smaller ones lie in parallel paths closer to the wall and farther upstream. Nose 
vortices are a particular example of secondary flow. With a round-nosed aerofoil, 
the boundary layer on the plate or wall and its vortex lines are rotated near the 
stagnation point, where the flow divides and turns through up to 90° in a very small 
distance. The secondary vortices due to this acute local curvature are very stable 
and persist for a considerable distance downstream. If the radius of the nose is 
reduced, then the position of the vortex moves nearer to the nose until, with a 
bi-convex or wedge aerofoil, the nose curvature is so small that the vortex is in- 
distinguishable with the scale of this experiment. Very close to the plate the 
boundary layer air does not reach the region of the nose at all. Smoke trails 
indicate that in this very slow moving air the filaments proceed in the nose direction 
to a point about 4 chord away when they suddenly turn through a right angle and 
fan out. 


The effect of the nose vortex on the cascade flow depends to a very large extent 
on the angle of incidence. Primarily, for small angles of incidence it results in a 
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Theoretical exploration of the effect of nose vortices on the spanwise velocity component. 
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vortex filament lying along the suction surface which opposes the main secondary 
vorticity and a vortex filament lying close to the concave blade surface which rotates 
in the same direction as the secondary vorticity. This is shown diagrammatically in 
Fig. 14, which also indicates the probable effect of the combined vorticity in creating 
an asymmetrical variation of spanwise velocity across the pitch. Such an 
asymmetrical pattern has been measured. 


4. Conclusion 


This investigation has shown that one of the effects of secondary flow within 
the passage between two compressor blades is to cause a local separation of the 
suction surface boundary layer adjacent to the cascade wall. The extent of the 
region affected cannot yet be estimated, but further experiments may enable an 
empirical approach to be used. As a result of this locally stalled area the cascade 
wall boundary layer is appreciably thickened, and farther downstream the deflection 
of the flow over the concave surface is reduced. Actuator plane methods of analysis 
also accentuate this effect, since they also predict a smaller deflection close to the 
wall. Estimation of the air angle leaving a blade row is hampered by the increase in 
thickness of the wall boundary layer and therefore the movement of the secondary 
flow pattern away from the wall. 


Local curvature of the streamlines near the nose of a blade or blunt body also 
generates a very stable secondary flow pattern, the vortices of which, while lying 
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parallel along the plate, pass right round the nose. They must therefore be 
regarded as superimposed on the main secondary flow vorticity. On the convex 
blade surface they oppose the secondary vorticity but are additive near the concave 
surface, in most blading configurations. Other opposite effects have been found at 
high negative incidences, but then local curvature in the opposite direction may 
predominate. 


All the features mentioned have been demonstrated experimentally, although 
they do not seem amenable to any analysis except a qualitative one. 
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. SUMMARY: Some new solutions for steady incompressible laminar boundary 

d layer flow, obtained by GéGrtler, have been used to test the accuracy of two 

methods which are commonly used to predict separation. A modification of 

Stratford’s criterion for separation is given in this paper and is probably the 

most accurate and the simplest of all methods at present in use. Modified 

numerical functions are also given for Thwaites’s method of predicting the main 

; ) characteristics of the boundary layer over the whole surface, which improve g 
4 the accuracy of the method. 


1. Introduction 


| Numerous methods are available for approximate predictions of laminar 
boundary layer separations. Of these, two, due to Thwaites”’ and Stratford™, are 

e fast becoming accepted as the most useful on the grounds of accuracy and simplicity. 
NOTATION 


x,z distances parallel and normal to the plane wall z=0 
u__ velocity component parallel to the wall 
U(x) the main stream velocity 


Xm the position at which U”(x,)=0, or xn=0 if U has not a true 
maximum 
U,, the main stream velocity maximum, i.e. U (xn) 
co 
8, the displacement thickness, | (1 
co 
6, the momentum thickness, [3 (1 - i) dz 


0 
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vy the kinematic viscosity 


L (m)=2 {[H (m) +2] m+1(m)} 
p_ the pressure of the fluid 
C, the pressure coefficient, =(p - p,)/(4eUn?)=1-—U?/U,,’ 


{(uw.) dx 


0 
A= Xo) dC, /dx] 
x, the exact position of separation 


X,,X,’ estimated positions of separation according to the methods of Thwaites 
and Stratford respectively 


X,,X,/ estimated positions of separation according to the modified methods 
A,B,C,K constants, appearing in various formulae 


C,, pressure coefficient at the estimated separation position of Stratford’s 
method 


Cy. pressure coefficient at the estimated separation position of the modified 
Stratford method 


s denotes values at separation 


w denotes values at the wall. 


2. Previous Methods 


Thwaites, considering the momentum equation only in the forms it takes at the 
boundary and when integrated across the boundary layer, obtained the relation 


Ud 
where m= -6, / ‘ (2) 
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Here U is the velocity at the edge of the boundary layer, v is the kinematic viscosity, 
8, and 8, are the displacement and momentum thicknesses of the boundary layer, 


5, (Ou 

L(m)=2 {m[H (m) +2]+/(m)}, . ‘ ; (5) 


and the suffix “‘w”’ refers to values at the wall, the normal to which is in the 
z-direction. 


By carefully examining the known exact and approximate solutions of the 
laminar boundary layer equations, Thwaites chose average values of the functions 
H(m) and I (m), and determined that separation should occur when m=0-082. He 
then found that L (m) is very nearly a linear function of m, and so chose 


L(m) ~ 0-45 + 6m. . (6) 
He then integrated (1), using (2) and (6), and found 


xz 


0 


From (2) and (7) it follows that separation occurs when 


dU 
= 
0-082 = - 6, ak 


Stratford divided the boundary layer physically into two distinct regions. In 
the outer layer, he showed that it was possible to consider separately the effects 
upon the velocity profile of pressure gradient in the absence of viscosity (as in 
inviscid theory), and of viscous forces in the absence of a pressure gradient (as in a 
Blasius layer). In the inner sublayer there is a balance between pressure gradient 
and viscous forces. Hence there is a transition across the layer, pressure forces 
being balanced by viscous forces at the wall and by inertia forces at the edge of the 
boundary layer. 


Mathematically, Stratford obtained conditions for separation which are 
asymptotically valid for sudden separation due to sharp pressure gradients. He 
then added terms of higher order, chosen empirically so as to fit known precise 
results which have large distances to separation. Thus, by interpolation, all possible 
values of separation distance were covered. His final result is that separation 
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occurs at a distance x which is a solution of the equation 


dC,\? dx? 140-144 
(x x)°C, 7°64 x 10-°(1+0°354) 2 140-46 icy (9) 
dx 
(x-x,) 


and U,, is the maximum velocity, which occurs either at the pressure minimum, 
X=Xm, or, if the pressure gradient is everywhere unfavourable, at the beginning of 
the layer, x=0. The distance x, allows for an initial favourable pressure gradient, 
and is defined as* 


Note that x,=0 when there is no initial favourable gradient. 


Stratford also gave a simpler, but less accurate, formula than (9), and this is 
in fact the solution of the boundary layer equations which is valid asymptotically 
for x—>x,, d?p/dx?=0 but dp/dx—>oo, namely 


(x-x)°C, (4) =7-64x10". . (13) 


The greatest difficulty facing these two authors must surely have been the lack 
of accurate solutions of the laminar boundary layer equations, particularly for 
separating boundary layers. For example, the only information available to 
Thwaites for the determination of the value of m at which separation takes place 
was as follows: — 


(i) The exact solution of Falkner and Skan“ for which the skin friction is 
everywhere zero has m=0-068. 


(ii) The exact solution for U=1-.x obtained by Howarth™ and confirmed 
by Hartree’ has m=0-084 at separation. 


(iii) For Pohlhausen’s® method, notoriously unreliable near to separation, 
m=0-157. 


(iv) Thwaites estimated that in Hartree’s™ calculations for Schubauer’s ellipse 
the value of m at separation “‘ would appear to be about 0-082.” 


From this inadequate information Thwaites chose the value m=0-082. 


*This notation is not quite the same as that used by Stratford. 
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Recently, however, many new solutions of the laminar boundary layer equations 
for separating boundary layers have been obtained by Gortler® by a new series 
method. By using only a few terms of the series, Gértler was able to obtain the 
distribution of skin friction along the surface for about 90 per cent of the way to 
separation, and he then used a finite difference method to complete the distribution 
and so obtained the separation point. Gdrtler’s method may thus be described as 
a rigorous approximation, in the sense that an exact solution could be obtained if 
it were carried out sufficiently accurately. 


It is important to note that Gortler’s method may lead in practice to a small 
systematic error in the predictions of separation. The curves of skin friction against 
distance along the surface bend down quickly towards the axis before separation, 
in the region where his finite difference procedure is applied. This corresponds to 
the singularity at separation (Goldstein), which implies an infinite slope in the 
skin friction curve. It seems that, because of this singularity, the skin friction as 
determined by an exact solution of the laminar boundary layer equations will tend 
to zero near separation more quickly than can be detected by the finite difference 
procedure used by Gortler. Now for the profile U=1-.x, for which an exact 
solution has been obtained, it is just in this way that Gortler’s result differs from 
the exact solution of Howarth. For x <0-100 the results are identical, but for 
x> 0-100 the skin friction tends to zero more rapidly in the Howarth solution, 
separation occurring at x =0-120 as against 0-126 obtained by Gortler. It is possible 
then that the results obtained by Gortler systematically overestimate the distance to 
separation by a small amount. This systematic error is too small to affect the 
value of Gértler’s work as a method for solving laminar boundary layers, 
but the fact that it is systematic does mean that his results must be treated with 
caution when used as data in trying to obtain a quick approximate method. 


This conclusion is further borne out by the fact that the discrepancies between 
Gortler’s results and those given by Thwaites’s method and by Stratford’s are 
much less for the two cases U=1- x? and U=1-.x* than they are for the other 
distributions considered by Goértler. These two cases have also been worked out 
by Tani®°), using the same continuation procedure as was used with great accuracy 
by Howarth for the case U=1-x. Gértler has obtained good agreement with 
Tani’s results; thus for these two cases there is a strong presumption that Gortler’s 
results are accurate. However, his other uncorroborated results may be subject to 
small systematic errors, which may well account for part of the discrepancy between 
his results and those of Thwaites and Stratford. 


The aim of the present authors has been to reconsider the methods of Thwaites 
and Stratford in the light of these new solutions. A new formula is obtained, 
similar to the simpler, less accurate of the two Stratford formulae (13), and this is 
found to be more accurate even than the full Stratford formula (9). However, since 
Stratford’s method can be used only for predicting the position of separation, it is 
suggested that when the boundary layer thicknesses or the detailed distribution of 
skin friction along the boundary are to be estimated a modified Thwaites’s method 
should be used. The essential modification is that separation is chosen to take 
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TABLE I 
SOLUTIONS FOR THWAITES’S METHOD 
1 | (1+x)7?} 0°159 0°147 0-091 0°158 
2 | (1 +.)-? 0:078 0-069 0-097 0:074 
3 0°223 0°209 9-091 0°221 
4 (1—x)? 0:067 0-061 0:094 0-065 
cos x 0°410 0°370 0°108 0°383 
6 0°401 0°363 0-120 0°372 
7 sin x 1-902 1-786 0-171 1-800 
8 1— x? 0-271 0°259 0-093 0:268 
9 0°462 0-440 0°104 0:449 
10 Ellipse 1:983 1-890 0:098 1-934 
1—x 0°120 0-116 0-086 0°123 


place at a slightly higher value of m, and hence small modifications are required in 
the prescribed functions /(m) and H (m). 


3. Thwaites’s Method 


Firstly, the accuracy of Thwaites’s method is tested in the light of the new 
solutions obtained by Gortler. Table I illustrates the results, in which the first 
seven velocity profiles listed are those considered by Gortler, the eighth and ninth 
are those considered by both Gortler and Tani, while the tenth and eleventh are 
respectively those of Hartree’s solution for Schubauer’s ellipse and Howarth’s 
linearly retarded layer. In the case of Schubauer’s ellipse the modified pressure 
distribution, for which Hartree found separation at x= 1-983, has been taken as the 
specified data. The values x, are the separation points as predicted by “‘ exact” 
theory, and x, are the predictions of Thwaites’s method. 


It is evident that Thwaites’s method systematically underestimates the distance 
to separation. For the seven Gortler results the root-mean-square error is about 
9 per cent, while for the other four results, which are known with more certainty 
to be accurate, it is about 4 per cent. The cause of part of this larger error may 
well lie in the possible systematic error in the Gortler results (discussed in Section 2), 
but part of it clearly lies in Thwaites’s method. If a larger value of m were chosen 
for separation, then the value of x at the predicted Thwaites separation point would 
be increased in all these cases. 


To get some idea of the increase required in the prescribed value for m at 
separation, it was decided to assume that equation (6) remained unaltered, and then 
to find the value of 


m, = 0°45 dU _ dx 
dx 
which would lead to the correct separation point for each of the velocity distributions 


considered. These values of m, are also given in Table I. As expected, they are 
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all greater than 0-082. At first sight it might seem reasonable to take the mean of 
these results, that is 0-105, but this would not be a fair method, for the results are 
not equally sensitive to changes in m,. For example, in spite of the fact that the 
predicted separation position is in error by only 6 per cent, the value of m, for 
the distribution U=sin x would need to be more than doubled in order to remedy it. 


After some thought it was decided that the best value of m, for all-round 
accuracy would be 0-090, and the new predicted separation points are given as 
x, in Table I. It should be noted that, although there is still a small systematic 
discrepancy for the seven G6rtler results, its root-mean-square value has been 
reduced from 9 per cent to 5 per cent, which is small enough to be ignored in 
view of the possible systematic error in GGrtler’s results. In addition the root- 
mean-square error for the other four results has been reduced from 4 per cent to 
2 per cent. A further small increase, say to 0-094, in the value of m, would further 
reduce the systematic discrepancy for the Gortler distributions to 3 per cent, but 
would increase the error in the prediction of the Howarth case to nearly 6 per cent. 
It was therefore decided that nothing was to be gained by this procedure. 


It remains to be asked whether this change in m, from 0-082 runs counter to 
any of the prescribed functions of Thwaites’s method. Clearly the function / (rm), 
which describes the distribution of skin friction,.and the function H(m), which 
gives the ratio of the displacement and momentum thicknesses, require modifications. 
Thwaites chose /(m) as a good average value from various known solutions, with 
1(0:082)=0. When m> 0-02 these known solutions leave a fair amount of scope 
in the choice of an average, and it is quite easy to choose an alternative average 
having / (0:090)=0. Similarly, the function H (m) needs to be defined for the slightly 
greater range of values of m, and in such a way that at separation it takes a 
prescribed average value. Thwaites chose the value 3-70, which is a weighted 
mean of the values 3:83 (Howarth), 3-19 (Hartree), 3-50 (Pohlhausen™) and 
4-03 (Falkner and Skan®’). 


TABLE II 
PRESCRIBED FUNCTIONS FOR MODIFIED THWAITES’S METHOD 
m l(m) H(m) m | H(m) 
—0:25 0°500 2:00 0-040 0-153 2°81 
—0:20 0:463 2:07 0-048 0-138 2:87 
—0:14 0°404 2°18 0-056 0°122 2:94 
—0°12 0°382 2°23 0-060 0°113 2:99 
—0°10 0°359 2:28 0-064 0°104 3°04 
— 0-080 0°333 2°34 0-068 0-095 3-09 
—0:064 0-313 2°39 0-072 0-085 3°15 
—0:048 0:291 2°44 0-076 0-072 3:22 
—0:032 0:268 2:49 0-080 0-056 3-30 
—0:016 0°244 2°55 0:084 0-038 3°39 
0 0:220 2°61 0-086 0:027 3°44 
0016 0°195 2°67 0-088 0-015 3°49 
0-032 0°168 0-090 0 3°55 
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TABLE Ill 
SOLUTION FOR STRATFORD’S METHOD 
1 0-159 0°150 0157 0-256 0-244 0253 
2 (1+x)-? 0-078 0-071 0-073 0-260 0:240 0245 
3 0-223 0-218 0-218 0223 0:218 0-218 
4 (1—x?? 0-067 0:063 0-064 0-242 0-229 0-232 
5 cos x 0-410 0-387 0-387 0-159 0142 0142 ?} 
6 1—x3 0-401 0-381 0-380 0-125 0-108 0107 
7 sin x 1-902 1-827 1-832 0-106 0:064 0-068 
8 1—x? 0-271 0-271 0-271 0-141 0-141 0-141 
9 0-462 0-462 0-461 0:0889  0-0889 00884 
10 Ellipse 1-983 1-976 2-01 0:0824 
1—x 0-120 0-120 0-121 0:226 0:226 0-227 


The present authors suggest that the value at separation of H=8,/5, should 
be determined in the following slightly modified way. In using Thwaites’s method, | 
5, is determined using (7) and (8), after which (3) is used to determine 4,. It 
follows that, if the method is to predict 5, as accurately as possible, it is best to 
choose H as the mean of ratios of exact 5, to approximate 5,. The values of H at } 
separation, as determined by this criterion, have been calculated using the new 
approximate values of 5,. They are 3°70 for the Howarth profile and 3-38 for 
Schubauer’s experimental distribution as considered by Hartree. Further, Gortler 
gives the displacement thickness for the velocity distribution U=(1+ x)~’, and his 


] 
value of 6, at sevaration would be predicted accurately by choosing H=3-:39. From ‘ 
these three values it would appear that H at separation should be rather less than ' 
3-70. After weighting in favour of the Howarth result the authors suggest that at , 
separation H is taken as 3-55. 

The suggested values for the functions /(m) and H(m), which do not require ; 
modification when m <0-06, are given in Table II. The modified approximate 
solution of the laminar boundary layer equations then becomes ; 

a 
dx a 
ou U fi 

4. Stratford’s Method ¥ 

As for Thwaites’s method, the accuracy of the method is tested against the n 
new solutions, and in Table III the results for the eleven velocity distributions under e 
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consideration are listed. x, is the “‘exact’’ separation point, and x,’ is the estimate 
—given by the full Stratford formula. Similarly C,, is the value at separation of the 
pressure coefficient C,, and C,, is Stratford’s estimate. 


Stratford’s method is rather more accurate than Thwaites’s method. For the 
first seven velocity distributions there is a systematic error in the predictions, as 
would be expected, but the remaining four cases are predicted with remarkable 
accuracy, both as regards separation point and the pressure coefficient at separation. 
In view of the fact that the root-mean-square value of the discrepancy in separation 
point for the seven Gortler profiles is less than 6 per cent, and Gortler’s values are 
subject to a possible systematic error of this order, little purpose will be served by 
trying to find a very much more accurate method than Stratford’s. The authors 
have therefore concentrated on the task of finding a simpler formula of comparable 
accuracy. In doing this they were prepared to accept a slight loss of accuracy in 
the last four cases, say 1 per cent accuracy (Stratford’s formula is more accurate 
than this in these four cases), while trying if possible to reduce slightly the systematic 
discrepancy in the first seven cases, say from 6 per cent to 4 per cent. 


Stratford’s simple formula, which is asymptotically valid for sufficiently sharp 
constant pressure gradients and sudden separation, is (neglecting the equivalent 
distance which must be included when there is an initially favourable gradient) 


(dG) _ a. 
#¢, (42) =164x10> 


By adding terms of higher order, and fitting them exactly to two known solutions, 
Stratford then obtained his much more complex formula (9), the additional terms 
representing the effects of varying distances to separation and varying pressure 
distributions. 


The present authors considered the possibility of interpolating, as Stratford 
did, with respect to A and d?C,/dx?, but in a much simpler way. Therefore they. 
considered an interpolation of the form 


d?C,/dx? 
(dC,/dx)* J 


2 
ve, (2) =7-64 x 10° {1+BA+C (20) 


dx 
and examined the values of B and C which would lead to exact answers for various 
velocity distributions. This type of formula is much simpler than Stratford’s (9), 
and is the simplest possible formula which includes both A and a non-dimensional 
form of d?C,/dx?. The task of choosing B and C so as to get good predictions of 
pressure coefficients, thereby automatically ensuring good predictions of separation 
points in all the cases considered, was very difficult. The method used was to 
work out the values of the various terms in equation (20) for each of the velocity 
distributions, thus obtaining a set of linear equations for B and C. These were 
then solved by the method of least squares, suitable weights being first given to the 
various distributions; for example, the last four distributions were weighted much 
more heavily than the first seven. It was found that the case U=sin x made these 
equations very difficult to deal with, and any attempt to improve even slightly the 
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accuracy of the estimate in this case led to a much greater loss of accuracy in the 
others. The reason for this is not difficult to find. For this velocity distribution, 
separation occurs very soon after the pressure minimum; hence, near to separation 
a small change in the value of x can lead to a large proportional change in the value 
of the pressure coefficient. It follows that the suspected error in Gortler’s results* 
will be greatly exaggerated in this case, and consequently his value of the pressure 
coefficient at separation is not to be relied upon. It was therefore decided to omit 
this case altogether, and B and C were chosen so as to get good all-round agreement 
with the other cases. The best values were found to be 


B=0°35, C=0-088 . ‘ (21) 


The predictions of the formula (20), with B and C given by (21), are of comparable 
accuracy to Stratford’s predictions and are more easily obtained. 


In carrying out this work it was found that the predicted separation points were 
remarkably insensitive to changes in B and C. For example, B and C can be 
changed by amounts of order 20 per cent and yet the predicted separation points 
vary by one per cent or thereabouts. It follows that in a formula of the form 


d’C, | (dC, 

| 

the results are more sensitive to changes in the value of A than to changes in the 

values of B and C. The authors therefore examined the possibility of obtaining a 

formula like (22) with B and C zero. The values of x*C, (dC,/dx)? were evaluated 

at the “‘exact”’ separation point for each of the cases under consideration, and an 

average value chosen, it being borne in mind (i) that weightings are to be given to 

the various cases on the grounds of accuracy and importance, and (ii) that the 

various cases are not equally sensitive to the value of A. The result of this choice 
led to 


@ 


2 
xc, (5) =104x10%. 83) 


Equation (23) can be easily solved for the various velocity distributions 
considered and the results are given in Table III. The values x,’ for separation 
are to be compared with x,’ (Stratford) and x, (“‘exact’’). Similarly the pressure 
coefficient at separation is given as C,., which is to be compared with C,, (Stratford) 
and C,, (““exact”’). In the prediction of the separation point the root-mean-square 
error for the seven Gortler values is 4-5 per cent, and for the other values it is 
0-8 per cent. Because of the possible systematic error of order 4 per cent in the 
Gortler values, the method appears to be exceedingly accurate, the error being no 
more than 1-2 per cent. For predictions of both separation point and pressure 
coefficient at separation, the method is at least as accurate as the full Stratford 
formula, is very considerably shorter, and is as simple in application as the simpler, 
less accurate of the two Stratford formulae. 


*Gortler found it more difficult to use his method for flows starting from a stagnation point, 
with a subsequent pressure minimum, than for flows everywhere retarded. 
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TABLE IV 
DATA FOR USE OF THE LIGHTHILL CRITERION 


Example Example 
1 (i+x)"! 0°86, 5 cos x 0-92. 
(1+ x)-? 0°86, 8 1-—x? 0:92, 
3 0-88, 6 0-94. 
11 (1— x) 0°87, 9 1—x* 0°95. 
4 (1—x)? 0:87, 10 | Ellipse 0-95, 


In view of the accuracy of the formula, at least in the pressure distributions 
considered, there seems no point in trying to refine the method any further by 
assuming non-zero values of B and C in (22). This would only lead to a more 
difficult formula to handle, and any gain in accuracy which might possibly be 
achieved would be quite unwanted. It is also doubtful whether any greater accuracy 
could be achieved by this means for an experimentally given pressure distribution, 
as it would require double differentiation of the pressure, as against the single 
differentiation in (23). 


The formula (23) and the simple Stratford formula (19) are both special cases 
of a more general formula which may be derived from considerations of the outer 
and inner parts of the boundary layer. It can be shown that 


{vo =x . 


where K is a function of the pressure distribution. Stratford’s simple formula has 
shown that, within the accuracy of his computations, K =0-00764 if the pressure is 
constant up to a position O and if, downstream of O, there is an asymptotically large 
constant pressure gradient. The work described here has shown that for less steep 
gradients of the type met with in practice a more accurate value of K is 0:01040. 
However, it is easy to construct pressure distributions for which K is very different 
from this, although such cases would be of academic rather than practical interest. 


5. A Note on a Criterion Applied by Lighthill 


In his paper on boundary layers and upstream influence, Lighthill“” considered 
a series of velocity profiles to which none of the usual methods can be applied, the 
main stream velocity gradient becoming infinite at separation. With a linearly 
retarded main stream, separation occurs when the main stream speed has fallen to 
0:88 of its initial value (Howarth), and Lighthill argues that the ratio 0-91 is more 
appropriate to the less favourable velocity curves that he considered. This raises an 
interesting point as to how the shape of the velocity distribution affects the fraction 
of the maximum velocity at which separation occurs. To help to settle this 
question, Table IV shows various values of U,/Um, as given by the present authors’ 
version of Stratford’s method. Here U, is the velocity at separation. 


For the first five cases listed, all of which are initially linearly retarded but 
differ as x increases, the flow separates when U,/U, is about 87-88 per cent. The 
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remaining cases, however, illustrate that the effect of a less favourable distribution 
is slightly more marked than Lighthill had suspected. For example, for U=1-— x’, 
where the gradient becomes progressively more adverse, separation occurs when the 
speed has dropped to about 93 per cent of its initial maximum value. For U=1-x', 
where the retardation is less in the early part and then later much greater, separation 
occurs before the speed has reached 95 per cent of its maximum. Thus it would 
seem at first sight that, for the profiles considered by Lighthill, 94 per cent would 
probably have been more realistic than 91 per cent. However, it must be 
remembered that Lighthill was interested not in the separation point precisely but 
in the point downstream of separation where the velocity gradient has fallen 
approximately to zero. With a real boundary layer the velocity here is rather less 
than at separation, because immediately downstream of separation the boundary 
layer can still withstand some pressure gradient. In using the Lighthill criterion it is 
necessary to include this effect in choosing the velocity fraction at separation; one 
must therefore be cautious in concluding from Table IV that Lighthill would have 
done well to have chosen a higher value. Certainly, if this criterion be applied 
in any future work, much more information regarding retarded boundary layers 
will then be available than was the case before Gortler’s results were published. 
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Maximum Ranges of Intercontinental Missiles 


DEREK F. LAWDEN, M.A. 


(Canterbury University College, New Zealand) 


SUMMARY : The thrust magnitude programme of a rocket missile being supposed 
specified, the problem is solved of programming the thrust direction to achieve 
maximum range on the Earth’s surface. Earth curvature and rotation, and 
variation of gravity with height are taken into account, but air resistance is 
neglected. Motion is assumed to take place in the Equatorial plane. Tables 
are included from which the optimal trajectories of missiles may be computed 
when two design parameters are known and these tables are used to assess the 
advantage to be gained at extreme ranges by firing in the sense of the Earth’s 
rotation rather than against it. 


1. Introduction 


In an earlier paper“, the problem of programming the thrust direction to 
achieve maximum range with a rocket missile has been investigated on the assump- 
tion that the effects of Earth curvature and rotation, air resistance and variation of 
the gravitational attraction with height above the Earth’s surface, could all be 
neglected. With an intercontinental missile, the air resistance will always be small 
by comparison with the thrust and weight forces and will rapidly become entirely 
negligible as the rocket ascends beyond a height of 20 miles or so. However, the 
effects of Earth curvature and rotation, and variation in g, must be taken into 
account or serious error will result. The object of this paper, therefore, is to. 
determine how the previous analysis must be modified to allow for these factors. 


NOTATION 
a=P/p 
A,B,C see equations (2) - (4) 
f component of missile’s acceleration due to the motor thrust 
g acceleration due to gravity at the Earth’s surface 
h=q,/V 


h (short-range missile) 


= h+w (long-range missile) 
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H height of the vertex of the missile’s trajectory 
k = g,/V 
l,,1, direction cosines of direction of motor thrust 
m_ mass ratio of missile 
M,_ mass of missile at launching 
motor thrust 
41:42 co-ordinates of missile at “ all-burnt ” 
G;.G2 Velocity components of missile at “ all-burnt ” 
r distance of missile from the Earth’s centre 
R radius of the Earth 
t time after launching 
T time of flight to “ all-burnt ” 


u=k+A/yv 

V circular velocity at Earth’s surface 

w=W/V 

W _ velocity of a point on the Equator due to the Earth’s rotation 
x=q,/R 


X,,X, co-ordinates of missile at time t 
X see equation (31) 
y=4q,/R 
a difference of longitude between point of impact and “ all-burnt” 
e angle of rotation of the Earth during time of flight of missile 
6 longitude of missile referred to a fixed meridian 

A,¥ see equations (25) or (29) 

# mass rate of consumption of propellant 
p angular range of the missile 
t= M,/p 
@ inclination of the rocket thrust to the horizontal 


2. Preliminary Results 


To ensure a two-dimensional problem, it is assumed that the missile moves in 
the Equatorial plane and we shall neglect the effects of Earth curvature and rotation 
and the variation in g during that phase of the motion when the rocket motor is 
operative. Taking horizontal and vertical axes Ox,, Ox., respectively in the plane 
of motion through the launching point O and rotating with the Earth, let (q,, ¢.) be 
the rocket’s co-ordinates and g,, q, be its velocity components at “ all-burnt.” 
Relative to fixed axes, instantaneously coincident with the positions of the axes 
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Ox,x, at “all-burnt,” the rocket’s velocity components parallel and perpendicular 
to the Earth’s surface are taken to be g, + W.q., where W is the speed of a point on 
the Equator due to the Earth’s rotation. If the rocket is fired against the Earth’s 
rotation, the sign of W must be reversed. 


At any instant during the subsequent motion of the rocket under the Earth’s 
attraction, let r be its distance from the Earth’s centre expressed in terms of the 
Earth’s radius R as unit, and let @ be the difference between its longitude at this 
instant and its longitude at “all-burnt.” The motion during this ballistic phase is 
supposed to be referred to non-rotating axes and longitudes are measured from a 
meridian plane fixed in space. Then, if V is the circular velocity at the Earth’s 
surface (i.e. the velocity of a body moving under the Earth’s attraction in vacuo in 
a circular orbit of radius equal to that of the Earth), the polar equation of the 
elliptical orbit into which the rocket is launched can be expressed in the form 


where, if y=q,/R, h=q,/V, k=q./V, w=W/V, then 


(3) 
k 


At the point where the orbit intersects the Earth’s surface, r=1 and @ satisfies the 
equation 


B+(1-—B) cos sin 0=A. ‘ ‘ (5) 
If « is the appropriate solution of this equation, then it can be calculated that 


C+{C?+(A-1)(2B-A-1)} 


and « is the difference in longitude between the missile’s point of impact and 
“all-burnt.” 


At the vertex of the rocket trajectory, it is found from equation (1) that the 
distance from the Earth’s centre is H, where 


A 
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The effective range of the rocket in angular measure can now be computed 
from « by subtracting the angular rotation of the Earth during the time « when the 
rocket is in free orbit and adding the increase in the rocket’s longitude relative to 
the Earth during the burning phase of its motion. The latter correction is 


(8) 


Since the moment of the rocket’s velocity about the Earth’s centre is constant during 
the ballistic phase, 


A (h+w)V 


R (9) 


r?6—constant = 


Hence, since on optimal trajectories r never takes values differing greatly from unity, 
and A=1 approximately, we have, very nearly, 


 (h+w)V 
A= (10) 
Thus the time taken to describe the ballistic part of the rocket track is given 
approximately by 


aR 


a) 


The angular velocity of rotation of the Earth being W/R, during this time < it will 
have rotated through an angle 


aw aw 


(12) 
The net angular range of the missile is accordingly 
w ah 


For a missile with a relatively small range, the approximation introduced into 
the calculation of ¢ will not be greatly in error. For a long-range missile, the error 
will be increased, but then the relevant term, namely w/(h+w), in equation (13) 
will be small, since w is then small by comparison with h. Hence, in either case, 
p will be subject to a small error only. 


The rocket thrust direction programme which maximises p is computed in the 
next section. 
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3. Optimal Thrust Direction Programme 


Experience suggests that the maximum range which is obtainable from a rocket 
missile is not critically dependent upon the thrust direction programme. For this 
reason, there is justification in adopting a relatively rough approximation for p 
before computing this programme. Thus, in equation (13), x is neglected altogether, 
since the horizontal distance travelled by the rocket during the operation of its 
motors will comprise a small fraction only (less than 10 per cent) of its total range. 
y, also, will be negligible by comparison with unity and hence equations (2) - (4) 
will be replaced by 


| k 


_ kth+w) 
Then equation (6) yields tan4a= (15) 


Also, for a long-range missile w/(h+w) is small and hence equation (13) may be 
written approximately 


« « « « 


and the problem has been reduced to that of maximising « (and therefore tan 4 a, 
as given by equation (15)). 


For a short-range missile, w/(h+w) will no longer be small, but the Earth’s 
rotation can now be neglected entirely, with only small effect upon the range and, 
therefore, even smaller effect on the optimal programme, and the whole motion can 
be referred to axes fixed in the Earth. Putting w=0 in equation (15), the angular 
range « from “ all-burnt ” relative to the Earth is given by 


hk 
tanda= (17) 
Also, equation (13) simplifies to p=2+ x, and thus, neglecting x, z, as given by the 


last equation, must be maximised. 


If /,, 1, are the direction cosines of the direction of thrust relative to the axes 
Ox,x,, it is shown in Ref. 1 that, provided that Earth curvature and rotation, 
g-variation and air resistance during the burning phase are neglected, to maximise 
the range p=p (q,, 2, G,, 72) the thrust direction programme must be chosen so that 


2 
(18) 

(T-)+5— 

0G, 


where T is the time of burning. 
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Hence, for a long-range rocket for which « as given by equation (15) is to be 
maximised, it is required that 


kfl+(h+w)’] 


For a short-range rocket, equation (17) is relevant and then 
k(+h?) 


Thus in either case the direction of thrust is constant. 


As in Ref. 1, we shall consider two cases, (i) that of a rocket whose thrust 
decreases exponentially in such a manner that the acceleration f it causes in the 
missile is constant, and (ii) that of a rocket whose thrust is constant. 


In case (i), as shown in Ref. 1, 
=41,fT’, q.=4(Lf—g)T. . & 


and hence, in terms of dimensionless quantities, 


1 

1 
where A=fT/V, v=f/zg 


and the relationship V?=gR has been used. fT is the velocity increment which 


would have been transmitted to the rocket by the motor thrust if this force were 


acting alone. A accordingly measures the fraction which this increment is of the 
circular velocity. If A=1 and the propellant is expended in negligible time, the 
rocket will achieve circular velocity. v is the acceleration caused by the motor 
thrust expressed in g’s. 


In case (ii), if P is the constant motor thrust, » is the mass rate of consumption 
of propellant and M, is the mass of the rocket at launching, the acceleration resulting 
from the thrust at a later time t¢ is f (t), where 


P a 
and a=P/p, r=M,/p. 
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Then, as proved in Ref. 1, the situation at “all-burnt” is determined by 
the equations 


q=al, [ log 
(28) 


2 
The dimensionless equations corresponding to equations (27) are identical with 
equations (23), provided that the parameters A and » are defined in this case by the 
equations 


a a T a 
A= log = Flog m, v= ef (29) 
where m=7/(7—T) is the mass ratio of the missile. The dimensionless equations 
corresponding to equations (28) are then 


1 1 


A table of values of X will be found in the Appendix. It will be observed that 
equations (30) reduce to equations (24) if X=4. Thus, with the appropriate defini- 
tions of the parameters A, v and X, equations (23) and (30) may be used in either 
case to determine the conditions at “ all-burnt.” 


Substituting for /, and /, in equations (23) the values of these direction cosines 
yielding maximum range (given by equation (19) for a long-range rocket and by 
equation (20) for a short-range rocket), results in the following conditions: — 


where h’=h for a short-range missile and h’=h+w for a long-range missile. A 
little manipulation shows that these conditions are equivalent to 


A 
k+ =V(Q?-h’) 
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If h’=h, it is easy to eliminate h between these equations and to show that 
u(=k-+A/y) satisfies the cubic 


ut+v (2 


Equations (34) and (35) have been solved numerically for and k for a selection 
of values of the parameters A and v. In each case /, and /, were computed from 
equations (23) and the angle ¢ made by the direction of thrust with the horizontal 
was determined. Tables of values of ¢ are given in Section 4. Given ¢, the extent 
of the rocket motion during the burning phase can be calculated from equations (30) 
and h, k can be re-calculated from equations (23). The range of the missile then 
follows from equation (13) and the height of the vertex from equation (7). 


4. Numerical Results 


Equations (32) and (33) were solved numerically in the cases when the design 
parameters took the values A=0, 0-2, 0-4, 0:6, 0:8, v=1, 3, 5, 7,9, 00. Two values 
of ¢ were computed in each case, one based upon the assumption that the missile was 
short-range and the other upon the assumption that it was long-range. These pairs 
of values were found to differ by less than 10 per cent in all cases and by less than 
one per cent at the intermediate ranges where the category to which the missile 
belonged was least obvious. No difficulty was experienced, therefore, in selecting 
an optimal value of # appropriate to each pair of values of A and v. This value was 
tested by computing the effect upon the range of positive and negative corrections of 
one degree. This effect was always found to be small (of the order of two minutes 
of arc) but, where necessary, adjustments were made to yield ¢-values, correct to one 
degree, as shown in Table I. Beneath each 9-value is shown the range along the 
Equator in degrees and minutes of arc. The value of A required to achieve circular 
velocity (A.) for each value of v considered is shown in the extreme right hand 
column. This was computed from the formula 


where w was taken equal to 0:0587 here as elsewhere. When computing the range, 
x and y were found from equations (24), rather than from equations (30). How- 
ever, the results will not be in error by more than a few minutes of arc in a case 
where equations (30) are appropriate. 

The calculations were then repeated with w= — 0:0587, i.e. assuming that the 
rocket was launched against the rotation of the Earth. The results have been 
collected together in Table II. 

It will be observed that at ranges greater than 2,000 miles (p > 29°), the 


difference between a missile’s range to the east and its range to the west is in excess 
of 10 per cent and, at ranges beyond 4,000 miles, this difference is more than 


25 per cent. 
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TABLE I 
OPTIMAL INCLINATION ~¢ OF ROCKET THRUST TO HORIZONTAL AND ANGULAR RANGE pp. 
(WITH THE EARTH’S ROTATION) 


A=0 0-2 0-4 os | 
90° 90° 90° 90° 
0° 0° 0° 0° 0 o | 
50° 49° 46° ti « 
4 49° 48° 4s° | 
f 48° 47° 45° 40° 
47° 46° 44° 40° 32° 
9 f | 6 
45° 44° 42° 37° 27° 
oo f 


TABLE II 


OPTIMAL INCLINATION ¢# OF ROCKET THRUST TO HORIZONTAL AND ANGULAR RANGE pp. 
(AGAINST THE EARTH’S ROTATION) 


A=0 0-2 0-4 0-6 08 10 | A, 
90° 90° 90° 90° 9° 
=" 1 oo 0 0° 0° oO 
5 f 49° 48° 46° 44° 39° 30° \ 1-081 
| 0° 2 8° 140? 42" 89° 19" f 
{48° 47° 45° 42° 38° 
7 oo 2° 4 88 38 ar |: 
0° 0 2° 8 52" 45° 54’ 101° 43" 
4s 45° 44° 39° 33° 
2° 207 9 40 23° aorta? tr | 10%? 


| 


As the velocity at “all-burnt” approaches circular velocity (i.e. at extreme 
ranges) the angle ¢ becomes small and the early part of the trajectory will pass 
obliquely through the atmospheric layer. Air resistance will then be an important 
factor in determining the shape of the initial arc of the optimal trajectory, causing 
it to rise steeply until the atmospheric film has been penetrated, after which its form 
will be governed by the equations of this paper. The calculation of the precise form 
of this initial arc will be the subject of further investigation. 


It will also be noted from Tables I and II that ranges greater than 4,000 miles 
correspond to a A-interval of length less than 0:2. Designing missiles to operate at 
all such ranges is an engineering problem of an order of magnitude comparable with 
that of the design of an artificial satellite vehicle. 


In cases where the optimal angle ¢ is small and the rocket thrust is constant, 
the reader is to be warned that the vertical component of rocket thrust may be 
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smaller than the initial weight of the missile. The early part of the trajectory will 
then intersect the Earth’s surface and the solution obtained is clearly quite 
unrealistic. This anomaly will disappear when air resistance is taken into account. 
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Appendix 
VALUES OF X GIVEN BY EQUuaATION (31) 
m 1:0 15 2:0 25 3:0 40 4:5 5:0 
X 0:5000 04662 0:4427 0°4246 0°4101 0:3982 0:3880 0°3792 0:3715 0°3645 
m 6°0 6°5 70 75 8-0 8:5 9-0 9°5 10:0 10°5 
X 0°3582 0:3524 0:3472 0:3425 0:3380 0°3340 0°3301 0°3265 0°3232 0:3199 
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Notes on the Mathematical Design of 
Nozzles and Wind Tunnel Contractions 


N. T. BLOOMER, B.Sc., Ph.D. 


(The University, Nottingham) 


SUMMARY: In this work the co-ordinates of three-dimensional nozzles of several 

types have been obtained as explicit functions, and some shapes have been plotted. 

It is hoped that some of these might be of use in cases where smooth uniform fluid 

flow is desired. The design of a particular nozzle shape in the presence of a plane 

boundary is given: previously any equations obtained have been too 
intractable to use. 


1. Introduction 


Some methods have been developed by Lilley“), Harrap(?) and others to obtain 
stream surfaces that have monotonic increasing (or decreasing) velocity (and hence 
pressure) distributions along the boundary walls and which would be suitable as 
nozzles or wind tunnel contractions. Hsue-Shen Tsien) has shown how velocity 
components “ and v and Stokes stream function ¥ may be obtained from a chosen 
axial velocity distribution, and use is made of this here to obtain the equations for 
nozzles of several types in an explicit form. His result is given as 


where f2"(x) denotes the 2n" differential coefficient of axial velocity f(x) with respect 
to x. 


An equivalent expression is 


u (x, So (x+ ir cos 6)d6. 
0 


1% 
Then and 


0 
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If f(x) is of the form A(x*+-a*)-" (where n is a positive integer), it may be written 
x [A, 
=0 


and repeatedly differentiated with ease. Further it was found that on substituting 
in (1) the resulting expansion is of binomial form and may therefore be summed. 


It should be clearly appreciated that stream surfaces along which the velocity is 
monotonic increasing or decreasing can only be obtained by trial and error. This 
process is not as laborious as it sounds, however, as any maxima in velocity dis- 
tributions are likely to occur near the singularities in the field. It is therefore advisable 
to plot the stream function near these regions first. 


NOTATION 
x co-ordinate measured along the axis of symmetry 
r co-ordinate measured normal to the axis of symmetry 
X,R_ co-ordinates of nozzles with unit throat radius 

¢ velocity potential function 
¥ Stokes stream function 

G, defined just after equation (5) 
G_ defined just after equation (7) 


pie parameters defined just after equation (17) 
constants 
=u along the axis of symmetry 
u velocity component in direction of x 
v velocity component in direction of r 
q_ resultant velocity 
Wall speed ratio 


axial speed ratio 


Suffixes 
w_ refers to velocities at the wall 


a_ refers to velocities along the central axis 
t refers to velocities at the throat 
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Stream function traces for certain configurations of singularities. 


2. Design of a Convergent-Divergent Nozzle 


By assuming an axial velocity distribution of 


it is possible to deduce the corresponding values of u, v and ¥, while suitable choice 
of b and ¥ enables one to plot stream surfaces with convenient wall speeds. This 
has been done. It is seen that a factor of the form (x?+1)-” in the axial velocity 
gives a ring singularity, of radius 1, in the plane x=0, as shown in Fig. 1. The axial 
velocity given by (2) will give surfaces as in Fig. 2. 


The equations deduced from (2), in the way briefly described in Appendix I, are 


2 


p=1 


3 (1 -G,*) sin yp 


+ 


3[x+(—1)?b]cos y,__ r? sin (3) 


2 
3siny, 3 [x+(—1)?b] cos3y, Ssin3y, , 3r* sin 5y, 
p=1 


2 


| 


3 
P G, 


3[x+( — 1)?d] siny,-3cosy, , 3r? cos 3y, 


+ Gs [x+-(— 1)?5]'sin Sy, cos Sr} (5) 
P 
P Pp 
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and G,t= + (— 


As an example, b has been set equal to 1/3, and —16)(1+5’)-* put equal to 
0-293278. 


A scale drawing and plots of the velocity ratios Q, and Q, against X appear in 
Fig. 3, together with some co-ordinate values of X and R. 


3. Design of Nozzles 


To obtain a nozzle of an asymptotic type, i.e. one in which the throat is at great 
distance, a configuration such as is shown in Fig. 4 may be used. As the singularities 
approach each other the axial velocity can be defined as 


mdb 
1e)= | 
0 


By setting n=2 and choosing m to give f,(00)=1, we have 


x 
10 
‘Be 
| 
0-75 6 3 
x 
@ VALUES 
0-2 
wx 
© 
% 5 4 3 r 1 0 
x 
Fig. 3. 


Showing the shape of a convergent-divergent nozzle and the magnitude of Q, and Q, at 
various values of X. 


x 0 OS 1:0 149 20 25 3005 35 40 425 450 4785 487 495 60 
R 1:00 1:00 1:02 1:04 1:10 1:18 1:29 1:52 1:81 201 2:26 268 290 346 — 
Q, 100 1:00 0:99 095 086 0:74 060 045 032 0:26 020 O15 O12 O10 — 
Q, 099 098 0:96 0:90 0:82 0:70 O57 045 0:34 0:29 O25 O21 0:20 019 0-095 
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AXIS OF SYMMETRY 


~ 


SINGULARITIES 


Fig. 4. 
Stream function traces for an infinite number of equal singularities. 


This is easily differentiated repeatedly and, upon substitution in (1), gives 
almost at once (see Appendix II), 


u (x, r)= {tan-4(G cos y) + (7) 
where 


G’sin 2y=2x, G’cos 2y=—1+7°+2’, and G4=4x?+(1 —x?—r*)?, as before. 


The value of —¥ is then found to be 


2 
(G c0s 9) - xGsiny G 
2 
and v= 1) sin y— x cos ‘ . (9) 


A scale drawing of the shape given by —’=0-1049 and plots of the velocity 
ratios Q, and Q, against X appear in Fig. 5, together with some co-ordinate values 
of X and R. 

4. Design of Wind Tunnel or Pipe Contractions 


Wind tunnel or pipe contractions can be obtained at once by the superposition 
of streaming parallel to the axis of symmetry. The axial speed may then be taken as 


and the constants adjusted according to the ratio required. 
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1-0 7 4 
|% 
le x 


© O,VALUES 
x 
q 


VELOCITY RATIO 
bd 


8 


Xx 
Fig. 5. 


Showing the shape of an asymptotic type nozzle and the magnitude of Q, and Q, at various 
values of X. 


co 6655) 64:37) «3:27 «218 «1°64 1:09 0°50 0 —0:38 -0°56 -—0:38 0 -1:09 
1:00 1:00 1:01 1:02 1:04 1:07 1:12 1:23) 141 1:72 2:18 2:66 2°82 


x 
R 
Q, 1:00 0:994 0:98 0:97 0:93 089 083 0:72 0:57 0-40 0°23 0-13 O11 
Q, 1:00 0993 0:98 0:96 0:91 086 0:78 0:64 0:50 0°40 0°34 0-40 0:50 0-23 0:00 


For a ratio of 3 : 1, say, 


2 x 2, 


A scale drawing of the shape given by —Y=0-2 and plots of the velocity ratios 
Q,, and Q, against X appear in Fig. 6, together with some co-ordinate values of 
X and R. 
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2 
0-75 Ox 
a™ 
z 0-50 x | 
fo) 
6 ° R 
0-25 
Q, VALUES 
© Q WALUES 
0} 
x 
Fig. 6. 
Showing the shape of the 3:1 contraction and the magnitude of Q, and Q, at various 
values of X. 
co 316 1:58 0395 O16 0 -016 -0-79 -1'58 -3:16 
R100 101 1:03 107) 217 1:22) 1-295 141 15400 
Q, 1:00 099 096 091 079 073 0-60 0°44 037 0340-33 
Q, 1:00 0:99 0:94 0:77 0:70 067 0-62 0°57 0-48 0-33 


a 


5. Equations relating to Nozzles in the Presence of Infinite Plane 
Boundaries 


The relative simplicity of the equations (7), (8) and (9) makes it possible to 
obtain expressions for a nozzle placed symmetrically in front of a plane boundary. 
By having two sets of singularities, of the same strength but different sense placed 
at a distance 2a apart with the same axis of symmetry, it follows that the flow through 
the symmetrical plane normal to this axis is zero. It may therefore be considered to 
be a plane boundary. 


The equations are then 


-1 -1 (x+ a) 
u (x, 0)= (x—a)+ tan (14) 
1 (r2 2 - a)G, siny, 
(G, cos y,)+ (G, cos y.) 


sin y, (G, (x- @) COS: sin y, (G, (17) 
1 2 
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(3) 
© ix 
© 
< ( 
= > 
0-503 2 
0257 
y 
fo Pe 
y, X Q,VALUES 
y VALUES 
1 2 4 6 7 8 


Fig. 7. 


Showing the shape of the nozzle in the presence of a plane boundary, and the magnitude of 
Q,, and Q, at various values of X. 


x co 9:03 682 5:72 462 407 352 2:93 242 205 1:92 205 242 165 1:10 0 
R 1:00 =1:00 1:01 1:02 1:05 1:07 1:13 1:24 1:44 1:77 2:20 2:55 2:74 
Q, 100 0:992 0:98 0:96 0:92 088 082 0:71 0:56 0:39 0:26 O18 O14 
Q, 100 0:99 0:98 0:95 0:90 084 077 0:63 0:50 038 035 038 O50 027 O15 0 


where G,’ cos 2y, = (x—a)*+r°—1, G,’ sin 2y, = 2(x—a), 
G;? cos (x+a)?+r?—1, G,? sin 2(x+a). 


The choice of a and the value of ¥ taken naturally decides the distance of the 
stream surface from this boundary. In the example given a=1-1 and —Y=0-103171. 
A scale drawing and plots of the velocity ratios Q, and Q, against X appear in 
Fig. 7, together with some co-ordinate values of X and R. 
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Appendix I 
CALCULATION OF u, v AND y GIVEN THAT {1+ 


It is only necessary to consider the case for fo(x)={1+ x*}-*, since the axial velocity 
specified can be obtained by the superposition of two similar terms 


i 3 


3i 


Differentiating 2” times with respect to x gives 


(QQn+2)!i  3(2n+1)! 


8[u"(x, 


+conjugate terms | 


Now H. S. Tsien’s result is 


u(x, r)= @) 


n=0 


which on substitution gives 


—1)"r"((2n4+2)!i 3(2 ! j ! . 
n=0 


This series is convergent for r?< |(x+i)?|, and on summation gives 


16u(x, r)= conjugate terms. 


By putting G? e*”=(x+i)?+,°, this gives 


from which equation (4) is readily obtained. 
The corresponding values of v and y may be found by application of 


0 
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Appendix II 
CALCULATION OF v AND yf GIVEN THAT tan i + 


This is considered in two parts: — 


(i) Consider first the case u(x, 0)=tan- x. 


1+x* 2i E i =¢ | 
1 1 1 


Putting these expressions into 


>< 1)" r"[u2"(x, 


2*"(n!)? 
n=0 
n=0 


the series being convergent for | (x + i)?|. 


It can be shown that the series 


1 (—1)"" r"(2n—1)! A+ +r?) 


and so, by putting x+ i=, x—i=a, u(x, r) can be expressed as 


2a 


Remembering that tan-! x5, log (3) +-K, say, 


where K is a constant, 


u(x, le + (i+ +K. 


B+ (6? +r’) 
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Setting a? +r =G? 


where G is defined as positive, substitution yields 


| 


cos y 


—tan- 
1+ Gsin y 


+C, say, 


where C is a constant. 


Since, however, when r is zero, u(x, 0)=tan-'!x, C is zero and, replacing x by $G?sin 2y, 
reduction gives 


u(x, r)=tan-! (G cos y). 
(ii) For u(x, 0)= 1 1 
rie xt+il’ 


On! 
so that (1) gives a(x, ‘ 
n=0 


which, on summation and the introduction of the parameters G and y, yields 


x 
It follows that, for fi(x)= [tn 1 


u(x, tan-(G 
(G cos y)+- G +3] 


The corresponding values of v and v v may be found by application of 


—Yy= and oy 
Ox 
0 


to the equation for u(x, r) and reducing as before. 
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Appendix III 


NoTE ON METHODS OF COMPUTATION AND POSITIONS OF SINGULARITIES 


It is possible to show that the positions of the traces left on a plane through the axis of 
symmetry by ring singularities in the field of flow, obtained by defining an axial velocity 
u(x, 0) are the same as the positions of singularities obtained in two-dimensional flow defined 
by the same axial velocity u(x, 0). 


The order of these ring singularities is rather less, however, than the corresponding 
two-dimensional ones, and so if the latter are very weak the former may not exist. It is 
nevertheless very instructive to consider the two-dimensional case so as to find the posi- 
tions of these singularities. The expressions obtained for u, v and y will then be valid 
for regions outside these ring singularities and rises in velocities along stream surfaces will 
generally occur near to any such singularities. Further they must clearly be excluded from 


the field of flow. 


The plotting of the stream surfaces is carried out by choosing axial values (say x = x;) 
and varying r until the correct value of y is obtained. Since in all cases the parameters for u 
and vy are the same, it is worthwhile each time to calculate values for u and apply the first 
approximation of 


u or. 


By straddling the value of ¥ required it is then possible to interpolate to whatever 
accuracy is required. 
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Approximation for the Effect of Twist 
on the Vibration of a Turbine Blade* 


A. I. MARTIN 
(Nottingham) 


SUMMARY: Modern turbine blades are generally twisted about their axes, the shape 
resembling that of a helicoid. This makes the determination of their overtone 
frequencies much more difficult. The equations governing the flexural vibrations 
may be obtained from the Clebsch-Kirchhoff theory for rods which are naturally 


twisted. By using these equations, this paper considers a first-order approximation 
for the effect of uniform twist. 


1. Introduction 


Two principal types of vibration of a turbine blade are flexure and torsion. If 
the blade is initially prismatic, the well-known Bernoulli-Euler theory for ordinary 
cantilever beams gives a simple formula for the frequency of any flexural mode and 
Saint Venant’s theory deals with torsion. If the blade has initial twist, which is quite 
usual in modern practice, pure torsion is hardly any different, but flexural overtone 


frequencies are changed appreciably, being increased or decreased according to the 
cross-sectional dimensions. 


Several authors have considered the problem of the vibration of a twisted canti- 
lever beam)-‘), Diprima and Handelman have obtained correct equations of 
motion and have proceeded to solve them by making use of the Rayleigh-Ritz principle. 
The fundamental frequency is thereby easily obtained, but for “‘ fairly high ” overtones 
the work would become far too laborious, since any approximating mode function 
for an overtone must be made orthogonal to the fundamental and any lower overtone. 
Troesch, Anliker and Ziegler have also obtained correct equations and, for an in- 
finitely thin cross section, the frequency equation is given in the form of an eight-by- 
eight determinant. Solutions of this determinantal equation were obtained by inter- 
polation on a sequence-controlled computer. Other authors have started out from 
an assumption such as a component of the bending moment about a principal axis 
being given in terms of the “‘ curvature” according to an equation of the Bernoulli- 
Euler theory for ordinary cantilever beams. The different methods of solution of 


* The work described here was done for the Bristol Aeroplane Company. 
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the final equations are greatly varied. Thus, White has considered the problem by 
means of an iteration process on integral equations and other writers in the field have 
indicated a solution through sets of simultaneous algebraic equations. Such procedures 
seem to require the use of some form of analogue computing machine to extract 
overtone frequencies, and the higher the overtone the more difficult does the task 
become. Experimental results have been obtained by such workers as Geiger and 
Rosard, and it has already been pointed out that there are certain breadth/thickness 
ratios for which the vibration becomes degenerate. 


The mathematical theory of naturally bent and twisted rods is given in Chapter 
18 of Love’s treatise’*’. For completeness and convenience we shall here first consider 
the determination of the equations representing the flexural vibrations of a uniformly 
twisted cantilever beam by using two facts of this theory. The equations may then 
be solved approximately by treating the twist as a first order perturbation. Let f be 
the ratio of the frequency of vibration of the twisted beam to the corresponding 
frequency of vibration of the untwisted beam, i.e. f is the “‘ correction factor” for 
twist. Then (for uncambered cross sections at least) fis, in general, an even function 
of the total twist ©, and a first approximate formula would be of the form 


f=H=1+K0’, . ‘ (1) 
where K is a function of the shape of the cross section. One of the objects of this 
paper is to obtain an expression for K. In particular, it will be seen that for a given 
overtone (or harmonic) K depends only on the ratio of the two principal moments 
of inertia of the cross section. If the effect of twist on the frequency of vibration is 
small over a workable range of values of ©, as is known from experiment for the 
fundamental, such a parabolic approximation (1) may be expected to fail except for 
exceptionally small values of ©. However, when the effect of twist is appreciable, 
as is generally the case for overtones, relation (1) may be expected to hold over a 
fairly wide range. 


NOTATION 
f ratio between frequencies of twisted and untwisted beams 
© total twist (in radians) 
K _ see equation (1) and Table I 
k,« components of curvature of strained central line of beam 


oyz “rotating” system of axes taken along principal axes of arbitrary cross 
section of beam 


OYZ axes fixed in space and parallel to principal axes of root section 
x distance along axis of beam as measured from fixed end 

6(x) (=kx), twist distribution (uniform) 
1 length of beam 


I,, 1, two principal moments of inertia of a cross section 
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n by E Young’s modulus 
have 
lures =dy/dx, etc. 


tract m_ mass per unit length 
task 
wad angular frequency 
ness 
pter see equations (11) 
‘mly t= (A,/r)' 
frequency numbers for corresponding untwisted beam 
f be 
ding (=1-875, 4-694, 7-855, 10-996, 14-137, . . . ) 
for p =tan}7 
tion 
yo, = cos sin (nr~), cosh (nr~), sinh (nr-*), respectively 
1 A, B, C, 
(1) F.GH see Table I 
this 
ven 
ents 2. Expressions for the Curvatures 
ie Let oyz be a system of axes which coincide with the principal axes of a cross 


for section, and let OYZ be axes fixed in space and parallel to the principal axes of the 
“ root section; so that 


Y=ycos¢# —zsing@ and Z=ysing +zcos@, . ‘ (2) 


if where 6=4(x) is the distribution of twist along the beam, x being the distance 
» measured from the root section along the undeformed central line. Referred to the 
axes oyz, the components of curvature « and «’ of the strained central line of the 
beam combine to form a vector lying along the binormal of the strained central line, 
the magnitude of this vector being equal to the curvature 1/» of this curve (see p. 382 
of Ref. 8). But, if Y and Z are co-ordinates of a point on the displaced central line, 
one of Frenet’s formulae for twisted curves* shows that the vector with co-ordinates 
ss |) dY/dx* and d*Z/dx* is of magnitude 1/p and lies in a direction at right angles to the 
binormal and the undeformed central line. Referring this vector to the axes oyz, it 
is found that 


K 


sin cos #, and cos sing; . (3) 


* See any textbook on differential geometry. To the present order of approximation the element of 
length along the deformed central line is 5x, and the tangent to this curve may be taken parallel 
to the undeformed central line. 
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or, in terms of y and z, differentiating the equations in (2) gives a 
_(d0\2 
dx® dx® “dx dx' “\dx 
ane he dt” dx dx? \dx 
These are the required expressions for the curvatures. When 6=0, they reduce to 
two equations given on page 396 of Ref. 8. 


3. Fundamental Equations of Vibration 

Consider a uniform twist distribution and put 6(x)=kx in (4), where k=O/I is 
a constant, / being the length of the beam. For pure flexure, the strain energy of the 
deformed beam per unit length of the undeformed beam may be taken to be* 


SED, 2+ x”, d 
where J, and J, are the two principal moments of inertia of a cross section and E is 


Young’s modulus. In terms of the co-ordinates y and z, the total strain energy is 


thus 
f 


+ 2ky'— dx+ iE (y"— — k*y)*dx, 
t 


where y’ = dy/dx, etc. When there is no taper of the cross sections, J, and /, are con- 
stants. At any instant of time ¢, the total kinetic energy of vibration is 


dy\? (dz\? 
+ 
0 
where m is the mass per unit length and is constant for homogeneous untapered 
beams. 


By Hamilton’s principlet the equations of motion, together with the boundary 
conditions at the free end of the beam, are obtained by forming the variation 
8 f(U—T) dt=0 subject to the fixed boundary conditions 


holding at the fixed end x=0 of the beam. We obtain, for the equations of motion, 
EI, (y"” — 2kz” — k*y") — kK*EI, (y" — 2kz' — k*y) — 
” 0? 
—2kEI, (2” + 2ky — Kz!) = —m 
} 


ands EI, (2"" + 2ky” — k*z")— + 2ky’ — k*#z)+ 
+ 2kKEI, (y" —2kz"—k*y')=—m | 


* Ref. 8, p. 395, equation (24). For pure flexure the torsional rigidity C may be assumed to be 
zero. 


t See, for example, Chapter 7 of Ref. 8. 
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and, for the free boundary conditions, 


— 2k2'()— =2"() + 2ky'O—K2()=0 


7 
and y" ()—2k2"()— k*y'D=2" 


To determine the natural angular frequencies p of a twisted beam, 0y/d7? and 
d*z/at? on the right hand side of (6) are replaced by —p*y and —p?z in the ordinary 
way, thus givinga pair of simultaneous ordinary differential equations, subject to boun- 
dary conditions given by (5) and (7). Since all the coefficients are constant terms an 
exact solution is possible, but it yields a complicated form of frequency equation (see 
Appendix I). This paper is concerned only with an approximate solution of the form 
y=yotk*y, and z=kz,, and valid for sufficiently small values of k in most cases. 


4. Existence of the Solution and Formal Equivalence of Twisted 
and Untwisted Beams 


Equations equivalent to (5), (6) and (7) may be obtained in a particularly compact 
form by applying the transformation (2). The resulting equations represent the 
vibration of a twisted beam referred to the fixed axes OYZ. Such equations may 
also be obtained, of course, by substituting (3), instead of (4), in the expression for 
the strain energy, and again applying Hamilton’s principle. 


Putting 6=kx in (2) and differentiating twice, 
Y"=(y" —2kz’ — k*y) cos kx —(z"+2ky’ — k*z) sin kx 
and = Z"=(y"—2kz'—k*y) sin kx+(z" + 2ky’ — kz) cos kx, 


and the method of direct substitution quickly gives the equivalent equations. If we 


use the notation of matrices and write u= Fal. then (5), (6) and (7) are transformed 


into 
d? 


u=u'=0 at x=0 and =0 at x=/, 


where the “‘ moment of inertia’ [(x) is given by the matrix 


I, 
74 
ty 


where I,=1f, cos* kx+ I, sin® kx, sin* kx+ I, cos? kx 


and = J,,=(/,— 4) sin kx cos kx. 
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1,, 1, and Jy, are the two moments of inertia and the product of inertia of a cross 
section about axes parallel to the fixed axes OYZ. McCann and MacNeal? have 
given an equivalent electrical circuit for the differential equations in (8). 


If (x) and u were scalars and not matrices, the equations in (8) would be identical 
with those representing the vibration of an ordinary untwisted cantilever beam having 
a variable moment of inertia (x) and a displacement u. Certain theorems holding for 
an untwisted beam may therefore be carried over to the case of a twisted beam. In 
particular :— 


For each value of I,, I, and k there are an infinite number of natural frequencies 
and corresponding solutions of the equations (5), (6) and (7), these solutions being 
sinusoidal varying quantities with respect to the time. 


There is, however, one theorem with which great care must be exercised if it is 
to be carried over; namely, that dealing with the question of degeneracy. For a given 
natural frequency, the scalar equation representing the vibration of an untwisted 
beam can have only one corresponding solution, whereas (8) may have two such 
independent solutions corresponding to the same frequency. For example, if /,= /,, 
the equations in (8) correspond to two identical sets of scalar equations and for each 
natural frequency there are exactly two independent solutions. If k=0, the equations 
(8) represent the two component flexural vibrations of an untwisted beam and there 
is a case of degeneracy if, and only if, the ratio /,//, is such that a frequency of one 
of the components is equal to a frequency of the other. 


5. Internal Resonance. Tendency of Certain Turbine Blades to 
Break No Matter What External Force is Applied 


Suppose that at each cross section of the beam a force acts in the direction of the 
principal axis oy with density P,=P, (x, t). On an actual turbine blade such a force 
may be supplied by a fluctuating gas pressure which is in excess of the temporal 
mean pressure and is caused by rotor blades passing stator blades. Then P,, must be 
added to the right hand side of the first equation in (6). Assume that in the ensuing 
motion the order of magnitude of the amplitude of z does not exceed that of y, and 
suppose further that the twist k of the beam is a small quantity. This will be valid in 
general, and the object is to show that such an assumption no longer holds in certain 
exceptional cases. Approximate equations of motion for the forced vibration of the 
beam are given by 


Oy 
El, y" +m Py (x,0) 
t 
and EI,z”""+m an —2kE(,+1)y” . . (10) 


These two equations represent two forced vibrations of an ordinary cantilever beam 
which moves under external forces of P,(x, t) and —2kE(/,+1,)y”, respectively. 
(Actually for (10) the boundary conditions show that there is also a bending moment 
equal to —2kEI,y'(/) acting at the end of the beam, but this makes no difference 
to the argument). 
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If P, is a sinusoidal varying quantity with respect to time and with frequency 
equal to a natural frequency of vibration of an ordinary cantilever beam of moment 
of inertia J, (and with mass per unit length m, etc.), then (9) shows that y ultimately 
becomes infinite, and ordinary resonance occurs. For a force P,, which contains no 
sub-harmonics of frequency equal to that of the ordinary cantilever beam, y remains 
finite but, generally speaking, it contains sub-harmonics having frequencies equal 
to those of the ordinary cantilever beam. (This will most certainly be the case if the 
force P, is suddenly applied and the beam starts its vibration from a state of rest"). 
Thus the forcing term on the right hand side of (10) will contain, in addition to the 
frequencies involved in P,, frequencies of the natural vibrations of the ordinary 
cantilever beam of moment of inertia /,. It follows that resonance of equation (10) 
can be of two types: either a natural frequency of vibration of the ordinary cantilever 
beam of moment of inertia /, is equal to a frequency in P, , or it is equal to a natural 
frequency of vibration of the ordinary cantilever beam of moment of inertia /,. 
Thus, if the slightly twisted beam is such that in the corresponding untwisted beam 
a“ broadside frequency ”’ is equal to an “* edgewise frequency,” then (10) will show the 
phenomenon of resonance, the corresponding solution z ultimately becoming infinite 
no matter what external force P, is applied. This is a kind of internal resonance. 


Returning to the exact equations of motion (6), then, according to Section 4, for 
each value of /,, J, and k there are corresponding solutions which are sinusoidal 
varying quantities with respect to time. Thus, just as in the ordinary way, the solution 
corresponding to an external force P, will always remain finite except for ordinary 
resonance. To bring this into agreement with the foregoing discussion, the following 
conclusion must result :— 


If any broadside force distribution is applied to give a forced vibration of a slightly 
twisted cantilever beam in which a broadside frequency is equal to an edgewise frequency 
for the natural vibrations of the corresponding untwisted beam, then the resulting 
edgewise amplitude is ‘* greatly in excess”’ of that of the broadside amplitude. 

Since a turbine blade is relatively thin, a vibration with “‘ large” edgewise 
displacements is expected to give rise to high stresses, and hence there will be a greater 
tendency for breakage whenever a broadside frequency is equal to an edgewise 
frequency in the corresponding untwisted blade. 


6. Approximate Solution of Equations (5), (6) and (7) 


Replacing 0’y/0t? and 0°z/0t° on the right hand side of (6) by —p*y and —p’z, 
an approximate solution is sought of the form 


y=y tk’y,, z=kz,, and A=Ay+k”A, 


where A= mp?/(E/,), and where yo, ¥;, Z;, Ag and A, are independent of k. If k=0, this 
solution reduces to that corresponding to an ordinary cantilever beam of moment 
of inertia J,. If a broadside frequency is equal to an edgewise frequency for the 
corresponding untwisted beam, Section 5 shows that a solution such as (11) is definitely 
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invalid when k0; otherwise (11) is valid for sufficiently small values of k by the 
general theorem of small perturbations. It should be observed that the broadside 
displacement y and the angular frequency p have been taken to be even functions 
of the twist kK, whereas the edgewise displacement z is an odd function of k, this being 
clearly required by physical conditions. 


The coefficients yo, y;, 2, Ag and A, are obtained by substitution and equating 
like powers of k. This gives the differential equations 


Yo" —ApVo= 9, 
— 1) 2" +2(2r+ 1) yo” 


wn Ns 
where r=/,/I,. Also, (5) and (7) give the boundary conditions 


Yor Ye =, at x=0, 


(13) 
and yo" =yo” =0, y,"=22,' + Yo. = — 39, = 2," =0, at x=/ 


The coefficient yp represents the mode of vibration of the ordinary untwisted 


cantilever beam. It is supposed that this is normalised, ie. [ye dx= 1. 


Integrating by parts four times and using the first equation in (12), together 
with some of the conditions (13), gives 


— Ags) dx = — +f (¥o"” — AoYo) ax 
0 0 


= —4y9 (I yo’ (D—2y0' (D) (I. 


Hence, on multiplying the second equation in (12) by yo and integrating, it follows 
that 


1 
A, = —2Ar+1) 2,” Yo Ax — 2(2r+ dx — (I) yo’ 
3 
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TABLE I* 

AF—BG+CH 
(T+), — 
2r 

C=(r+IQE+ (Z—o), H=E+ rt 


Solving the third equation in (12) subject to the corresponding boundary con- 
ditions in (13) and substituting the result in (14), the correction factor f is obtained 
in the form f2=A/A,=1+k?A,/A,. Thus, to the present order of approximation, 
K=,/(2PA,) in (1). The work is straightforward but rather long. It is first observed 
that —2[(r+1)/(r—1)](1/A.) y»” is a particular integral of the equation for z,; it 
then remains to obtain the complementary function, say V(x). Writing y=cos tl, 
o=sin tl, ¢t] and where t=(A,/7)'. it is found that 


where (1+ Vix)= — r+) cos tx-+ 


Yo" (0) 


By using the two identities y?-++-7? = 1 and [?—X?=1, it is easy to verify that z, is 
indeed the required solution of the differential equation satisfying the correct boundary 
conditions. 


The condition 1+ y['=0 means that there is a broadside frequency equal to an 
edgewise frequency in the corresponding untwisted beam. If r+1, this means that 
the function z, becomes infinite, and the state of affairs described in Section 5 exists. 


When evaluating the integrals in (14), use has been made of entry (3) in Table II 
of Ref. (11), together with results in Appendix II of this paper. The values of y, and 
its derivatives at either end of the beam are given by equations (14) in Ref. 11. 


* For even harmonics; for odd harmonics replace p by —1/p, change the signs of the first terms in 
A, B and C, and the signs of rt in F and G. 
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Writing "—A,'/ for the frequency numberst of the corresponding untwisted 
beam and p=tan}", so that y=cos(nr-*), o=sin("r-*), and 
xX =sinh (r-*), the final result is given in Table I. 

At a first glance it might appear that z, and K(r) become infinite as r approaches 


unity. In fact, as r approaches unity, K(r) approaches zero, and z, approaches the 
limiting function / [4 — 2/(»,2p?) — x// ]yo. For this, the reader is referred to Appendix LII. 


As r->oo, y and [' approach unity, and * and & approach zero as fast as r=. 
Thus, A=2r+ O(r'), with a similar expression for B. Also, F=r'!+ O(1), with a 
similar expression for G. Clearly, C= O(r*?) and H= O(r-*). Hence 

AF —BG=[2r+ O(r')][r! + O(1)] — [2r+ O()][r! + O(1)] = O(r); 
and CH= O(r'). 


Thus, AF—BG+-CH=O(r), and the last term in the expression for K tends 
to. zero. Hence 


x) ~ S210) 


as r—>oo. For odd harmonics replace p by —1/p in this formula. A list of approximate 
numerical values for K is given in Table III of Section 8. 


7. Direction of Oscillation of the Tip of the Blade 


Let ® be the angle between the direction of oscillation of the tip section and the 
perpendicular to the root section, measured positive when taken in the same direction 
of rotation as the twist. Then, by (11) and (15), 


2). 
Thus, to the present order of approximation, 
Vil) 
o= 
and the appropriate values of V(/) and y,(/) are substituted. For even harmonics it is 
found that 


tan(®—@)= 


(r—1)(1+yI) np 


where A, B and C are the same expressions as in Table I. The corresponding formula 
for odd harmonics may be written down from the rule given in the footnote of Table I. 


t For the first few frequency numbers of an ordinary cantilever beam see Table I of Ref. 11 or any 
textbook on vibration theory. For the nth overtone (or nth harmonic). (n > 1) we have n=(n+4)= 
approximately. 
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For the fundamental (= 1-875) and r=4, ®/0 was calculated to be about 0:25. 
For the second overtone (77-855) and r= 100, ®/© was about 1-35. These values 
are in fair agreement with experimental results obtained by Geiger" using indicator 
card traces. 


8. Numerical Results 


A broadside frequency is equal to an edgewise frequency whenever r=(1/"')', 
where 1) and 1’ are any two frequency numbers for the ordinary cantilever beam. 
If r takes one of these critical values, the solution in Table I is no longer valid*, 
and the corresponding slightly twisted blade is more likely to fail in the engine. If 
r lies near to one of these critical values, then formula (1) will only hold for exception- 
ally small values of ©. 


The calculations are shown in Tables IJ and III. 


For high overtones and not too large values of r(e.g. 3rd overtoneand4<r< 144), 
land & become sensibly equal, and it then happens that BG and CH are of the same 
order of magnitude. To compute their difference it is then preferable to use a formula 
such as 


2r 


— yto i 


= B+ | (r+ 1)—2/r (y H, (odd harmonics), 


both sides being equal to B(G— H)+(B—C) H if we take [= >. 


TABLE II 


CRITICAL VALUES OF r(>1) FOR WHICH THE SOLUTION IN TABLE I FAILS 


Ist Overtone 
or 
1st Harmonic 


2nd Overtone 


or 
2nd Harmonic 


3rd Overtone 
or 
3rd Harmonic 


4th Overtone 
or 
4th Harmonic 


* In fact, according to general perturbation theory the common frequency will be split by the twist 
into two component frequencies lying on either side of it, the amount of splitting being pro- 
portional to the twist © and not to ©? as in (1). 
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TABLE Ill 
APPROXIMATE VALUES OF K= K (r) 


Ist Overtone 2nd Overtone 3rd Overtone 4th Overtone 
r or or or or 
Ist Harmonic 2nd Harmonic 3rd Harmonic 4th Harmonic 
1 0 0 0 0 
+0-084 +3-58 —0-86 —0-063 
9 +0-14 —9-86 —0-059 —0-028 
16 +0-254 —1-11 —0-044 —0-023 
25 +0°753 —0-461 —0-036 —0-02 
36 +3-14 —0-264 —0-04 —0-0167 
49 — 1:34 —0-174 —0-036 —0-0123 
64 —0°66 —0-126 —0-034 —0-0008 
81 —0-47 —0-097 —0-033 +0-312 
100 —0-38 —0-076 —0-032 —0-041 
121 —0-33 —0-06 —0-031 —0-035 
144 —0-046 —0-03 —0-03 
169 —0:274 —0-032 —0-0284 —0-0282 
196 —0-26 —0-015 —0-0272 —0-0268 
225 —0°244 +0-011 —0-0263 —0-026 
256 —0-233 +0-064 —0-0255 —0:0252 
289 —0-224 +0-303 —0-0251 —0-0246 
324 —0-217 —0-502 —0:0246 —0-0242 
361 —0-213 —0:198 —0-024 —0-0238 
400 —0-206 —0:135 —0-0237 —0-0234 
625 —0-19 —0-088 —0-0144 —0-0221 
1,296 —0-17 —0-075 —0:17 —0:0195 
2,401 —0:16 —0-07 —0-051 —0-0112 
4,096 —0-154 —0-068 —0-043 —0-037 
6,561 —0°151 —0-066 —0:04 —0-027 
10,000 —0-148 —0-065 —0-038 —0-025 
—0-137 —0-06 —0-034 —0-022 


A line between two “ frequencies ’’ in a column denotes passing through a critical value of r (Table II); 
K then changes sign from + to — 
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Appendix I 


EXacT FoRM OF THE FREQUENCY EQUATION 


The main paper has shown how to obtain a first approximate value for an overtone 
frequency of a twisted cantilever beam. To obtain a more exact value such an approximation 
may be substituted in a part of the exact form of the frequency equation and a solution 
obtained for the resulting expression (method of iteration). 


Replacing 0*y/dt? and 0*z/df? in (6) by —p*y and —p’z, it is required to solve the resulting 
pair of differential equations subject to conditions (5) and (7). It is simplest to work in terms 
of a Laplace transformation (see Ref. 12). We seek for solutions of the form 


where the integrations are taken round a simple closed curve lying near infinity in the com- 
plex €-plane, and where »(£) and {(£) are analytic functions regular and vanishing at £=00; 
i.e. these functions have Laurent expansions of the forms 


valid for all large values of £. By (5), 
[nce de =fen(e ae=0, 


with similar equations involving {(¢). Thus a,=a,=b,=6,=0; i.e. and have zeros 
of the third order at £=0. 


Substituting (16) in the pair of differential equations, we obtain 
and — + 21,)€? — (u — k*1,)} + 2k, + ef dE =0, 


where »=mp?*/E. The two expressions in the square brackets have simple poles at =o. 
Also, they cannot contain a series of descending powers of £, since such a series would make 
the corresponding integral non-zero. It follows that these expressions are linear functions 
of €; i.e. 


{1,£4 — + 21, — (u — k41,)} + bE 
and — 2k°(1, + 21,)€? —(u — K*1,)} + 2k, + LL =c + dé 


where a, b, c and d are constants. 


(17) 
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Putting (16) in the boundary conditions (7), we have 


and — + 2hEn( = KANE) + =O. 


(18) 


Solving (17) for »(€) and () and substituting the result in (18) gives four homogeneous 
equations for the determination of the constants a, b, c and d. For a non-zero solution the 
corresponding determinant of the fourth order must then vanish, and we obtain the frequency 
equation for the determination of the numbers pn. 


Thus, writing 


+ 4k°(, + 


the solution of (17) gives 
A, — k?)n(§) — 2KE(E)} = (a + BENE? — {1 (E? + - + + LT + + 
and 


A, — + (a + (EF + k?)? + w} + (c+ dEME? — (EF 


Let An= 5 and (n=0, 1, 2, 3, 4), 


so that A, and B, are real quantities (u real) satisfying An=Bn44+2k*Bn42+k*Bn. Then the 
frequency equation can be set in the form 


| 1,(A’,—k* A’) — k*B’») 2k(hA’,+ (19) 


=0, 


where A’» and B’m are two sub-matrices of the second order given by the equations 


, Am, Am+1 Bm, Bn +1 
A'm= , and B’n= > (m=0, 1, 2). 
Am +1, Am+2 Bm +1, Bm+2 


Since the number 4 is contained in every A and B it would be a difficult matter to solve 
(19) in the general case. However, knowing an approximate value of u, this equation may 
be used to obtain a more exact value by placing the approximate value in the relations 
defining A, and B, and solving the resulting equation for p. 


If k=O, the left hand side of (19) reduces to the product of two determinants of the 
second order. Placing these two determinants equal to zero separately, and evaluating the 
integrals defining A, and B, according to the calculus of residues, we obtain the well known 
frequency equations for an untwisted cantilever beam corresponding to moments of inertia 
I, and J,, respectively. 
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TABLE IV 
Integral Even Harmonics Odd Harmonics 
— 2 2)» 
(st—t) 
2 2),—1 2 2 
cosh tx dx - s(T—s*)p s(t? I'+s de} 
(st—t*)/1 
| (x) sinh tx dx 2{ =e} 


0 
s=n/l, p=tan}n, y=costl, o=sin tl, T'=cosh tl, X=sinh tl. 
y(x)=mode of vibration of snsitied cantilever beam of length / when the frequency number is 1. 


Appendix II 
SOME FURTHER INTEGRALS RELATING TO THE VIBRATION OF A CANTILEVER BEAM 


According to (14) and (15), it is necessary to evaluate the integral i (x) yo dx. Thus 


we require integrals of the modes of vibration of an ordinary euntiuees beam with the 
elementary circular and hyperbolic functions. The notation of Part I of Ref. 11 will be used. 
Thus, if y=y(x), 0«<x~«</, represents a mode of vibration, then y=K,(x)—K.(x), where 
K,"(x)= —s*K,(x), K,"(x)=s*K,(x), and s=»/l. It is easy to evaluate integrals of K, and K, 
with circular and hyperbolic functions. 


1 
Thus (s?—1*) K,(x) cos tx dx= cos” tx—K,"(x) cos tx}dx 
0 0 


{K,(x) cos’ tx—Ky'(x) cos tx}dx 
0 


= (x) cos’ tx —K,’(x) cos tx] 
0 


and the right hand side may be evaluated by formulae (9) in Ref. 11. There are seven 
similar expressions for the integrals of K.(x)cos rx, Ki(x) sin tx, K.(x) sin tx, K,(x) cosh tx, 
K,(x) cosh tx, K,(x) sinh tx and K,(x) sinh tx. Integrals such as 


| y(x) cos tx dx 


are then obtained by subtracting appropriate expressions. In this way Table IV may be 
constructed. 
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Appendix III 
THE CASE WHEN r APPROACHES UNITY. VIBRATION OF A SLIGHTLY TWISTED SQUARE PRISM 


The usual method of obtaining the limiting form, as r approaches unity, of an expression 
similar to that given in Table I is to expand the various quantities in powers of r—1. However 
in the present case the work becomes much too involved, and it is better to proceed in a 
slightly different way by first showing that the function z, approaches a limit which is then 
substituted in (14). 


Placing the values of A, (=n*/I*), yo” (0) and yo'(l) [see formulae (13) in 
Ref. 11] in the equation for V(x), gives 


(20) 


for even harmonics. At r=1, the denominator of this expression has a zero of the second 
order, since 1++I then reduces to the frequency expression for the ordinary cantilever 
beam. To deal with the numerator, use the identity 


which may easily be verified by multiplying out the expression on the left hand side and 
using the equations y*++-o?=I'*—=X?=1. By formulae (6) in Ref. 11, the second factor on 
the left hand side does not vanish at r=1 (even harmonics) and it follows that the first 
factor has a zero of the second order at this point. Hence, z,”(0) approaches a finite limit 
asr—>l. 


Similarly from the identity 


(r+ 
8roX(i+yP), 
we may obtain 
C+ oT) +2Vr(o—%) 
plirt (r—1) 


V” (0) 


and, by again making use of formulae (6) in Ref. 11, it easily follows that V” (0) is asymptotic 
to —(8y?//#)/(r—1), as r—>1. But, by (15), 


r+1 


21" )=-2 


yo"(0)+V" (0), 


and (14) of Ref. 11 showsthat the first term on the right hand side is asymptotic to (8y*//8)/(r— 1). 
Hence, by expanding in powers of (r—1), z,:”(0) is seen to approach a finite limit as r—>1. 


Since the boundary values z,(0)(=0), z,’(0)(=0), z:"(0) and z,”(0) are sufficient to 
determine completely the solution of the third equation in (12), and since this solution is a 
continuous function of these boundary values and of r, it follows that z,(x), together with 
its derivatives, tends to a finite limit as r—>1, the limits of the derivatives being the derivatives 
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of the limit. Hence we may make r—>1 in the third equation in (12) and in the corresponding 
boundary conditions in (13) to obtain 


with boundary conditions z,(0)=z,’(0)=0, z,”(D= 2y,’(I), z,’’"()=0, where in these 
equations z, now represents the value of this function at r= 1. 


It is easily seen that the most general solution of the differential system (21) is given 


where a is an arbitrary constant. 
Making r—>1 in (14), we have 


0 0 


where z, is now given by (22). It is seen that the terms corresponding to ay, cancel. Thus 
these terms contribute 


1 
—4a Ye" Yo dx—2ayy'(I)*=4a 
0 0 


since yo (0)=0. All the other terms also cancel. Thus their contribution is 


xyo" ye dx-+6 yo"yo dx—2ya(l) y+ 
0 0 


where, to obtain the second line, entries (3) and (4) of Table II and equations (14) of Ref.11 
have been used. It thus follows that K(r) approaches zero as r approaches one. Thus there is 
no first-order correction for the frequencies of a slightly twisted square prism. 


To obtain the limit of z, as r—>1 it is necessary to evaluate the constant a. This may be 
done by obtaining the actual value of the limit of the expression on the right hand side of 
(20). Alternatively, if we wish to avoid all actual expansions of quantities in powers of 
(r—1), then we proceed to obtain an integral condition. 

Integrating by parts, 

1 
Ag ” , ” Ay 
0 0 


and in this equation we substitute the values of y,”” and z,”” as given by the first and third 
equations in (12), together with appropriate boundary conditions in (13). Thus 


2/ 1 yo ax= *) yodx. 
0 0 
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The integral on the left hand side is given by entry (4) in Table II of Ref. 11; and substituting 
the values of yo(/) and X, gives 


/ 


for even harmonics. 


Thus 


provided that r-41, and by passing to the limit it follows that this integral condition also 
holds when r= 1. 


Putting (22) in (23) and evaluating the integrals (entries (1) and (11) in Table II of Ref. 11), 
we obtain 
a=(3— 


for even harmonics. For odd harmonics replace p by —1/p. 


If we proceed to determine K(1) in a direct manner without passing to a limit, then we 
are confronted with an indeterminate differential system (21) which has no unique solution 
unless some additional restriction is imposed, such as (23). However, as before, all solutions 
of (21) lead to K(1)=0. Thus, K(r) is continuous at r=1, while the corresponding edgewise 
displacement z,=z,(x, r) is discontinuous, having two linearly independent values at this 
point. This is, of course, typical of degenerate vibrations. 
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